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Os| ■ ABSTRACT 

| The effects of wind-driven star formation feedback on the spatio-temporal organi- 

I- zation of stars and gas in galaxies is studied using two-dimensional intermediate- 

representational quasi- hydrodynamical simulations. The model retains only a reduced 
subset of the physics, including mass and momentum conservation, fully nonlinear 
fluid advection, inelastic macroscopic interactions, threshold star formation, and mo- 
mentum forcing by winds from young star clusters on the surrounding gas. Expanding 
shells of swept-up gas evolve through the action of fluid advection to form a "tur- 
qq ■ bulent" network of interacting shell fragments whose overall appearance is a web of 

qq | filaments (in two dimensions). A new star cluster is formed whenever the column 

^C) . density through a filament exceeds a critical threshold based on the gravitational in- 

■ stability criterion for an expanding shell, which then generates a new expanding shell 

| after some time delay. A filament-finding algorithm is developed to locate the potential 

ON ■ sites of new star formation. 

ON ■ The major result is the dominance of multiple interactions between advectively- 

distorted shells in controlling the gas and star morphology, gas velocity distribution 
fS and mass spectrum of high mass density peaks, and the global star formation his- 

tory. The gas morphology strongly resembles the model envisioned by Norman & Silk 
(1980), and observations of gas in the LMC and in local molecular clouds. The depen- 
dence of the frequency distribution of present-to-past average global star formation 
rate on a number of parameters is investigated. Bursts of star formation only occur 
when the time-averaged star formation rate per unit area is low, or the system is small. 
Percolation does not play a role. The broad distribution observed in late-type galaxies 
k>( \ can be understood as a result of either small size or small metallicity, resulting in 

j_j ■ larger shell column densities required for gravitational instability. The star formation 

' rate increases with density, but dependences on gas velocity dispersion and average 

shell column density suggest that the dependence is multivariate. The distribution of 
gas velocities exhibits exponential tails over a broad range of parameter values and 
the velocity distribution for gas in filaments is nearly exponential. Decay simulations 
with no star formation suggest that the exponential tails are due to multiple shell 
interactions, not individual stellar winds. The cloud mass spectra, estimated using a 
simplified version of the structure tree method, tend to be power laws at the higher- 
mass end, with an index that is nearly independent of the star formation activity or 
model parameters. Kinetic energy decay in simulations without star formation yields a 
t^ 1 dependence. We discuss how the simulations can be viewed in the context of vari- 
ous incomplete conceptual models, including collisional cloud coalescence, wind-driven 
turbulence, propagating star formation, forced mass-conserving Burgers turbulence, 
and granular fluids. 

Key words: hydrodynamics, turbulence, stars: formation, ISM: bubbles, galaxies: 
ISM 
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1 INTRODUCTION 

The collective behavior of star formation in disk galaxies 
appears to involve interactions between many stellar and 
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interstellar processes acting over a large range of tempo- 
ral and spatial scales. Central to understanding the star 
formation process is the question of how spatial and ve- 
locity structure is generated in the interstellar gas through 
the action of instabilities, magnetic fields, turbulence, star- 
formation, shocks, etc., how that structure leads to the cre- 
ation of dense star-forming gas clouds, and how the resulting 
young stars act as feedback to the gas structure. The present 
paper attempts to address this question through the use of 
an "intermediate-level representation" of the hydrodynam- 
ics, in which gas structure and star formation are controlled 
by the interactions of wind-driven shells subject to nonlin- 
ear advection and a prescription for star formation in shells. 
These interactions lead to a "turbulent" evolving network 
or web of filaments (in two dimensions). 

There is now considerable observational evidence sup- 
porting the view that some form of "turbulence" is an impor- 
tant organizing process at work in the ISM over a wide range 
of spatial scales, although the physical nature and source of 
this "turbulence" remain unclear. Molecular line widths are 
usually much broader than the predicted thermal widths, 
indicating that the gas motions are supersonic (e.g. Larson 
1981). Furthermore, the line widths scale with the size of the 
region as a power-law (Larson 1981; see Elmegreen & Fal- 
garone 1996 for more current results), reminiscent of incom- 
pressible turbulence or a field of discontinuities, although 
striking exceptions exist (e.g. Plume et al. 1997). Evidence 
for a turbulent ISM is found at very small scales in pul- 
sar observations (Armstrong et al. 1995), and turbulence is 
often argued to provide a source of pressure in supporting 
the galactic disk (e.g. Ruzmaikin et al. 1988, McKee 1990, 
Lockman & Gehman 1991, and Falgarone et al. 1992). The 
observed gas morphology is also suggestive of a dynamic 
process at work. Molecular line, IRAS, HI, and extinction 
maps reveal irregular, often filamentary, scale-free gas dis- 
tributions on scales from 0.01 pc to hundreds of parsecs (e.g. 
Scalo 1985, Falgarone & Phillips 1991, Elmegreen & Falgar- 
one 1996, Chappell & Scalo 2000). This complex appearance 
suggests that the gas is in a highly agitated state and has 
not relaxed into equilibrium structures. A common interpre- 
tation of this observation is that the gas is "turbulent." 

A likely organizing process at work in the ISM is the 
action of high-mass stars and protostellar winds on the sur- 
rounding gas and on future generations of star formation. 
Winds and supernovae lead to the formation of shock waves 
which expand into the ambient gas. If the shocked gas has 
time to cool, a compressed gas shell forms behind the shock, 
which sweeps up the interstellar gas as it expands outward. 
Such shells may drive gas motions and turbulence in the ISM 
if their kinetic energy and momentum is efficiently trans- 
ferred to the surrounding gas. For example, Ruzmaikin et al. 
(1988) showed that the power input from supernovae may be 
sufficient to account for the overall interstellar cloud veloc- 
ity dispersion. Abbot (1982) and Castor (1993) estimated 
and compared the energy input from various massive star 
sources in the Milky Way, and Tarrab (1983) did the same 
for nearby galaxies. Norman & Ferrara (1997) calculated 
in more detail the spectrum of energy input from H II re- 
gions, supernovae, and superbubbles. This idea is the basis 
of several self-regulating star formation scenarios. Norman 
& Silk (1980) suggested that, for low-mass star forming re- 
gions, this energy could maintain turbulence in the parent 



cloud, preventing an overall gravitational collapse and lim- 
iting star formation. They viewed the turbulence as con- 
sisting of a network of old shell fragments described by a 
cloud fluid and argued that star formation may persist if 
the shell fragments coalesce to form more massive gravita- 
tionally unstable clouds, similar in spirit to the larger-scale 
simulation models to be presented here. Observational sup- 
port for this view is found in regions of vigorous star for- 
mation, in which the gas is often found to be organized 
into complex, possibly turbulent, networks of H II region- 
driven, wind-blown, or supernova shells (e.g. Braunsfurth & 
Feitzinger 1983, Meaburn 1984, Bruhweiler et al. 1991, Chu 
& Kennicutt 1994, and Kim et al. 1998 for studies of the the 
Large Magellanic Cloud; Chini et al. 1997, Bally et al. 1999 
for studies of the Orion and Circinus clouds, respectively). 

Expanding gas shells may also mediate propagating star 
formation. According to Elmegreen and Lada (1977), the 
shells may become gravitationaly unstable and lead to the 
formation of a new stellar generation. (See Vishniac 1994, 
Elmegreen 1994, 1999, and references therein for a more 
detailed discussion of the problem.) The winds and super- 
novae from the new stars propel the process such that star 
formation may propagate through a cloud or galaxy. Their 
model was originally developed to describe the spread of 
high-mass star formation through a molecular cloud. Prop- 
agating star formation theories have since been developed 
for shells driven by protostellar winds from low-mass stars 
(Norman & Silk 1980) and for large-scale propagating star 
formation mediated by superbubbles created around entire 
OB clusters (e.g. Elmegreen 1994, Jungwiert & Palous 1994, 
Silich et al. 1996, Ehlerova et al. 1997). It has also been 
suggested that star formation could proceed at the intersec- 
tions of converging flows in the turbulence as suggested by 
Elmegreen (1993) in a slightly different context. These sce- 
narios suggest a positive feedback mechanism in which star 
formation leads to the birth of more stars. Observational 
evidence has been found for propagating star formation as- 
sociated with superbubbles (Oey & Massey 1995, Lortet and 
Testor 1988, and Dopita et al. 1985), shells created by indi- 
vidual supernovae (Elmegreen 1989), and the compression 
of clouds on small scales by strong UV radiation fields from 
massive stars (Sugitani et al. 1995). Comprehensive surveys 
of the evidence are given by Elmegreen (1992, 1999). 

These observations suggest a picture of the ISM in 
which shells of swept-up (i.e. compressed) ambient gas cre- 
ated around sites of vigorous star formation interact to form 
a turbulent background of old shell remnants. This picture 
has observational support in regions of active star forma- 
tion. For example, in the giant HII region 30 Dor in the 
Large Magellanic Cloud, Chu & Kennicutt (1994) describe 
the morphology as a "complex network of expanding sys- 
tems" and cataloged many shells over a broad range of spa- 
tial scales and velocities. The H I aperture synthesis map of 
the entire LMC by Kim et al. (1998) is dominated by fila- 
ments, shells, and holes. The detailed study of the smaller- 
scale Circinus molecular cloud by Bally et al. (1999; see also 
Dobashi, Sato, & Mizuno 1998) illustrates how the same 
type of dynamical morphology occurs due to "churning" by 
outflows from lower-mass protostars. In addition, shell in- 
teractions may be the primary mechanism controlling star 
formation in these regions. For example, Bruhweiler et al. 
(1991), and Chu & Kennicutt (1994) report evidence that 
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star formation occurs at the intersections of supershells in 
30 Dor. Furthermore, Hyland et al. (1992) suggest that star 
formation is maintained in this region because of shell inter- 
actions. 

The recent HST/NICMOS study of 30 Dor in the LMC 
by Walborn et al. (1999), in conjuction with the morpho- 
logical evidence in Chu & Kennicutt (1994) and Kim et al. 
(1998), strongly suggests that shell networks and propagat- 
ing star formation are intimately related. This relation is 
one focus of the present work. This scenario of star forma- 
tion mediated by a network of shells was first proposed by 
Norman & Silk (1980) for low-mass star-forming regions, 
and the observational results of Bally et al. (1999) for the 
Circinus molecular cloud strongly resembles this picture. 
While Norman & Silk (1980) incorporated star formation 
feedback processes, they modeled the turbulence through a 
one-zone cloud coagulation model which neglects all non- 
linear spatial hydrodynamic interactions and assumes equi- 
librium solutions which rely on balancing source and sink 
terms in a set of ordinary differential equations. Such "one- 
zone" models are unable to investigate either spatial or dy- 
namic phenomena. The inclusion of nonlinearities, time de- 
lays, or multiple interacting "phases" into one-zone models 
frequently results in nonequilibrium behavior such as oscil- 
lations, bursts, or chaos. Nonequilibrium models have been 
applied to cloud fluid systems (Scalo & Struck-Marcell 1987, 
Vazquez & Scalo 1989), propagating star formation modeled 
as diffusion (Shore 1981, Ferrini & Marchesoni 1984, Kor- 
chagin et al. 1988), supernova-driven multicomponent ISM 
models (Ikeuchi & Tomita 1983, Ikeuchi et al. 1984), and 
self-regulating star formation systems (Parravano 1989). An 
important question regarding these models is the extent to 
which the resulting nonequilibrium behavior depends on the 
lack of spatial degrees of freedom in the model equations. 
More recently, Norman & Ferrara (1996) demonstrated the 
importance of star formation feedback by calculating the 
source spectrum of forcing due to HII regions, supernovae, 
and superbubbles. However, in the treatment of spatial cou- 
pling, they modeled transfer between scales in Fourier space 
by the advection term (we use "advection term" to refer 
to the u • Vu operator in the fluid momentum conserva- 
tion equation) through a transfer function applicable for 
incompressible flows, an assumption which contradicts the 
observed supersonic motions in the ISM. Furthermore, the 
spatial distribution was only treated at the level of the two- 
point correlation function. 

The main goal of the present work is to use simulations 
to study the spatial behavior caused by highly compressible 
advection in the presence of threshold star formation and 
stellar wind-driven expanding shells, and to investigate the 
global behavior that results from "opening up" the spatial 
degrees of freedom. We want to understand how stellar out- 
flows "churn" (Bally et al. 1999) the ISM. 

Several past studies have attempted to investigate the 
interaction between star formation and gas structure using 
hydrodynamic simulations. Chiang & Prendergast (1985) in- 
troduced a simulation in which stars are also modeled as a 
fluid, with the star formation rate proportional to the gas 
density to some power. In their simulations, a network of 
gas filaments surrounding regions of active star formation 
emerges as a consequence of a star-gas instability. Bregman 
et al. (1993) considered a similar model but also investigated 



a simulation oriented vertically in the galactic plane with 
the inclusion of a galactic gravitational potential. Vazquez- 
Semadeni et al. (1994) and Passot et al. (1995) studied a 
turbulent model which includes the effects of self-gravity, 
magnetic fields, rotation, stellar heating, and cooling and 
which models star formation as a threshold function of the 
gas density. Navarro & White (1993) also adopted a thresh- 
old star formation law in galaxy formation models. They 
found that their results depended sensitively on the pre- 
scription chosen for stellar energy injection. 

While full hydrodynamic approaches provide the most 
detailed information on the star-gas interaction, large simu- 
lations including many physical processes require supercom- 
puters and even then only a modest number of simulations 
can be run, and for a relatively short time. In complex prob- 
lems that contain many underdetermined parameters, large 
regions of parameter space must be left unexplored. Further- 
more, even if we were confident in the selection of parameter 
values and we somehow knew that the hydrodynamic sim- 
ulations accurately reproduced the true star-gas dynamics, 
we would still not be guaranteed a satisfactory theory of star 
formation in the sense of an explanatory model. As Kaneko 
(1991) puts it "If one succeeds in reproducing the phenom- 
ena from the model equation, then what can one learn? One 
dangerous trap in computational physics is that one may 
be still in the same level as the direct observation of the 
complex phenomena itself. What one might obtain is just 
that the equations are correct and reasonable, without any 
understanding of the complex phenomena." In astrophysics 
the situation is a bit different since we do not have the lux- 
ury of directly observing the "phenomena" over time. Thus, 
simply being able to watch the simulated evolution of the 
ISM, for example, could in itself be very useful. However, nu- 
merous examples from fields including non-linear dynamics, 
turbulence, biology, ecology, and economics, have found that 
in studying complex systems which exhibit self-organization 
or the emergence of non-trivial correlations, the most de- 
tailed simulations do not always provide the greatest insight 
into the phenomena. These studies find that experimenta- 
tion with the level of representation of a simulation can lead 
to new insight into the nature of the underlying dynamics 
and can play an important role in formulating new theo- 
ries. A recurrent finding of these investigations is the result 
that the phenomenology, statistical behavior, scaling laws, 
etc., of many "complex" systems can be reproduced by sim- 
ple models which possess the symmetries, conservation laws, 
nonlinearities, or spatial connectivities of the physical sys- 
tem without reference to the detailed form of the system. 

One of the more successful examples of this approach 
was used by Yanagita and Kaneko (1995) to study Rayleigh- 
Bernard convection. They combined a simple donor-cell ad- 
vection scheme with "rules" to simulate the effects of buoy- 
ancy and incompressibility and reproduced nearly all the 
phenomena observed in Rayleigh-Bernard convection exper- 
iments including the formation of convection rolls, oscilla- 
tions, routes to chaos, and spatio-temporal intermittency. 
In fact, since their method is computationally efficient, they 
were able to explore the parameter space of the model, and 
discovered a new class of phenomena. The present paper is 
similar to this convection study with respect to the level of 
representation. 

A number of papers have tried to model the spatial as- 
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pect of the galactic star formation problem by circumvent- 
ing the hydrodynamics, particularly in the area of propagat- 
ing star formation. Perhaps the best-known is the stochastic 
self-propagating star formation (SSPSF) model (see Seiden 
& Gerola 1982 for a review) in which active star-forming 
sites spread stochastically over a spatial lattice. While this 
model abstracts the detailed physics of the propagation 
mechanism to a simple stochastic rule, it proposes a pic- 
ture of the organization mechanism through the theoretical 
framework of percolation theory. Gerola and Seiden used 
this model to make predictions about the coherence of star 
formation in disk galaxies, in addition to global star for- 
mation histories in dwarf galaxies. Their model has been 
extended to include anisotropic effects (e.g. Jungwiert & 
Palous 1994) and 3-D geometries (Comins 1983). Lacking 
from the SSPSF picture is the potentially important role 
that the gas dynamics (in particular the nonlinearity of the 
advection operator and consequences of conservation laws) 
may have on the stellar organization. Also lacking from 
such models is an explicit treatment of the physics of the 
propagation process, which is simply imposed as a stochas- 
tic "rule." More ambitious cellular automaton-type mod- 
els have been presented by Perdang & Lejeune (1996) and 
Lejeune & Perdang (1996), who have generalized the Gerola 
and Seiden approach to include a highly deformable grid 
in the presence of a prescribed velocity field with both ro- 
tational and non-rotational (diffusive) components, as well 
as several stellar evolutionary states. Particularly interest- 
ing is the ability of these models to produce spatial struc- 
tures with fractal dimensions similar to those observed in 
the interstellar medium. However these models are still non- 
hydrodynamic and do not allow for any feedback between 
star formation and the (prescribed) velocity field. An even 
more detailed cellular automaton model for propagating star 
formation was presented by Gardiner et al. (1998), who in- 
clude an N-body calculation of the gravitational interac- 
tions and star formation feedback. However the simulations 
are still non-hydrodynamic, and so the important effects of 
nonlinear fluid advection are omitted. Another recent non- 
hydrodynamic cellular automaton model, concentrating on 
radiative transfer coupling between cells, was presented by 
Rousseau, Chate, & Le Bourlot (1998). 

Reaction-diffusion models of propagating star formation 
have also been studied. Shore (1983) introduced a model 
similar to the Lotka-Volterra equation in which the propa- 
gation of star formation is assumed to be controlled by a 
diffusion process. Feitzinger (1984), Neukirch & Feitzinger 
(1988), Neukirch & Hesse (1993), Nozakura & Ikeuchi (1984, 
1988), and Korchagin & Ryabtsev (1992) have used simi- 
lar reaction-diffusion descriptions to study the problem of 
pattern formation in the ISM. Patterns formed in reaction- 
diffusion models are often referred to as "dissipative struc- 
tures," following Nicolis & Prigogine (1977). The main draw- 
back of these models is their neglect of the underlying gas 
dynamics, in particular the replacement of the nonlinear ad- 
vective spatial operator with a linear diffusion operator for 
the gas. This severe assumption suppresses whole classes of 
phenomena (e.g. shocks) and may inappropriately introduce 
others (e.g. Turing pattern formation). For this reason we 
believe that recent re-introductions of this model (Smolin 
1996, Freund 1997) should be regarded with caution. 

In this paper, we develop a model for stellar "churn- 



ing" of the ISM which is structurally simpler than a full 
hydrodynamic simulation, but which retains the fundamen- 
tal mass and momentum conservation laws and the nonlinear 
hydrodynamical effects of fluid advection. However, pressure 
is neglected, magnetic fields are ignored, and the effects of 
self-gravity are abstracted as a threshold rule for star for- 
mation based on a linear gravitational instability criterion. 
Instead of attempting to build a "complete" model of the 
ISM, we have chosen to investigate this hybrid or "reduced" 
model to study the organization and evolution of star for- 
mation driven by the interactions of a system of turbulent 
wind-blown shells. While the structure of this model is in 
a sense simplistic, it includes enough physics to be of inter- 
est from several theoretical perspectives including theories 
of propagating star formation, turbulence, self-regulation, 
cloud coalescence, and pattern formation in general. 

The model developed in this paper simulates a system 
of interacting shells driven by winds and supernovae from 
massive stars, although the general results should be ap- 
plicable to the case of clouds whose internal structure and 
motions are driven by protostellar or protocluster winds. 
Stars are born from the shells when the gas column density 
through a shell exceeds a threshold value, and their winds 
generate new expanding shells, perpetuating the propaga- 
tion of star formation. The shell dynamics are controlled by 
fluid advection and obey global mass and momentum con- 
servation. The simplest set of fluid equations satisfying these 
requirements are the continuity equation and the pressure- 
less momentum equation. By neglecting pressure gradients, 
collisions between gas structures become completely inelas- 
tic, causing the colliding shells to "stick." This limit cor- 
responds to a vanishing adiabatic index 7 which has been 
estimated to be considerably less than unity in the inter- 
stellar medium, at least at low and moderate densities (see 
Vazquez-Semadeni et al. 1996 and Scalo et al. 1998) . It can 
also be interpreted as assuming that ram pressure dominates 
thermal pressure everywhere, i.e. that the motions are every- 
where highly supersonic. (However these two interpretations 
are not equivalent, as recently shown by Passot & Vazquez- 
Semadeni 1998.) By neglecting the pressure gradient term 
and by only incorporating self-gravity in terms of a threshold 
condition on star formation based on a gravitational insta- 
bility condition, the computational speed increases consid- 
erably, allowing us to better explore parameter space and 
test variations of the model. In addition, the simplification 
of the model equations allows for easier interpretation and 
opens up possibilities for constructing analytic theories of 
the dynamics. 

We refer to this model as a "reaction-advection" sys- 
tem since the gas motions are driven by young star clusters 
(the formation of which depends on the gas properties) and 
evolve through the action of nonlinear fluid advection. We 
suggest that, like reaction-diffusion models, it may be gen- 
eralized to study generic spatio-temporal pattern formation, 
but would be applicable to systems in which the transport 
process is best described by fluid advection rather than dif- 
fusion. The model is, in effect, equivalent to forced Burg- 
ers turbulence, although our "forcing" is tied to the physics 
of star formation and stellar winds rather than simply as- 
suming some arbitrary (usually Gaussian) stochastic forcing, 
and we explicitly solve the mass continuity equation. 

In section 2, we present the fluid equations, star forma- 
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tion criteria, and stellar forcing prescription, and discuss the 
computational methods employed. In section 3, the spatial 
distribution is presented for model runs which include star 
formation feedback and for decay runs with no stellar forc- 
ing. Section 4 discusses the global properties of the simula- 
tions, including the distributions of present to past average 
SFR ratios, the gas velocity, and the cloud mass spectra. 
An interesting aspect of the simulations is that they can 
be viewed in terms of several more "highly-reduced" mod- 
els such as collisional cloud coalescence or propagating star 
formation. A discussion of this "multiperspectival" aspect is 
given in section 5, followed by a summary. 



2 THE MODEL 

The motions of the interstellar gas in our model are con- 
trolled by fluid advection and momentum forcing due to star 
formation: 



dp 



at 

dpv 



+ V ■ (pv) = 



p w N*(x',t) 



(1) 
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where p is the gas surface density, p m is the total momentum 
input per massive star, iV*(x',i) is the number of stars per 
unit area injecting momentum at position x', and t w is the 
duration of the momentum injection (10 7 yr here). Star for- 
mation is assumed to occur according to a shell column den- 
sity threshold criterion meant to represent shell gravitational 
instability, as explained in section 2.2 below. In practice we 
also experimented with models that include the source term 
in the continuity equation (1), accounting for the gas lost to 
star formation. However because the timescale for this gas 
depletion is so large compared to the phenomena of inter- 
est, and in order to calculate averages and other statistical 
quantitiesfrom a stationary distribution, we have omitted 
the mass depletion term for the calculations reported here. 

The momentum source term in eq. 2 only symbolically 
represents the finite difference procedure followed in the sim- 
ulations. The position in question is x. The sum is over the 
eight nearest neighbor cells, at positions x'. The unit vector 
ensures that the momentum is directed toward the position 
x. If there is a cluster at x', then Ni.(x') is the number 
of newly-formed stars at that position, per unit area. The 
number of stars formed in that cluster is computed from the 
mass in the cluster using an adopted IMF. The cluster mass 
is computed from the mass in the simulation cell times an 
assumed constant star formation efficiency. The momentum 
input p w = m(x')v is calculated using a constant assumed 
wind velocity, 40 km s _1 , at a distance corresponding to a 
cell size (7.8 pc in the simulations reported here), and m(x') 
is the fraction of the mass released by the cluster at x' that 
enters the cell at x, assuming the morphology of the wind is 
circular in two dimensions. (We do not yet address the in- 
teresting question of the effect of collimated outflows rather 
than spherical winds.) The division by t w signifies that this 
mass and momentum are redistributed over the lifetime of 
the wind, 10 7 yr. Momentum is only injected after a time 
Td since the onset of gravitational instability (see section 2.2 
below). The motion of the cluster at position x' is taken into 



account when calculating velocities. Star formation is in ef- 
fect externally specified by selecting a criterion, generally a 
function of p and v, as discussed below in section 2.2. 

The effects of pressure gradients were excluded from the 
model for both computational and conceptual reasons. The 
relative importance of pressure to advection in the ISM may 
be estimated through dimensional analysis as 



/' 



'VP 7c 2 L" 



7 

M 2 



v-Vv " v 2 L~^ " M 2 ^ 
where c is the adiabatic sound speed, L is the characteristic 
scale of the flow, and M is the characteristic Mach number. 
Since, as discussed above, 7 may be small and M is large 
in much of the ISM, the effects of pressure are frequently 
dominated by advection. Several observational studies have 
found that macroscopic "turbulent" gas motions in the ISM 
dominate thermal pressure, especially at large scales (e.g. 
Falgarone et al. 1992) and that gas motions in the interstel- 
lar medium are frequently supersonic with Mach numbers 
from a few to 100 (e.g. Larson 1981) so the neglect of ther- 
mal pressure may be a good approximation. The utility of 
the Burgers equation for studying highly supersonic astro- 
physical flows has been demonstrated by Kofman & Raga 
(1992 and references therein) in the context of stellar jet 
structure. Still, it should be remembered that, with constant 
pressure, the system has some unique properties: e.g. there 
are no sound (pressure) waves to carry information and all 
compressions develop into discontinuities due to the lack of 
pressure gradient forces. Passot & Vazquez- Semadeni (1998) 
have used a model equation to show that the zero-pressure 
behavior is not equivalent to the limit 7^0, although the 
major difference they find is for the lowest- density regions, 
which are not of primary interest in the present work. 

In addition, the effects of pressure cannot be neglected 
near shock discontinuities. In the present simulations the 
widths of filamentary structures formed by shocks and their 
interactions are instead controlled by the (unphysical) nu- 
merical diffusion, which sets the filament widths initially 
to a few cell widths and causes them to gradually fatten 
with time (see sec. 2.1 and Appendix A). Such an effect 
could mimic the effects of real turbulence within the fila- 
ments (if it exists), but the scaling of this physical effect 
is unknown, and so our filament thicknesses cannot be con- 
sidered realistic. On the other hand, the critical quantity 
for the star formation threshold that we adopt is the col- 
umn density through a filament at a given location, and this 
quantity should not depend much on the filament thickness. 
The identification of filaments and the calculation of their 
local column densities is described in sec. 2.3 below. 



2.1 Donor-cell advection 

We choose the donor-cell scheme to model fluid advection 
on a discrete lattice. The main advantage from our perspec- 
tive is its straightforward representation which may be easily 
adapted to include momentum injection by star formation 
and, in a future study, moving grids to simulate galactic dif- 
ferential rotation. The donor-cell method may allow higher 
dynamic ranges in the fluid density than, for example, SPH 
codes (Durisen et al. 1986). However, the use of adaptive 
smoothing lengths in modern SPH codes reduces this dis- 
parity (see Kang et al. 1994). Also, since the method is 
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explicit, expensive matrix inversion techniques are not re- 
quired as in implicit methods (van Leer 1977). Donor-cell 
methods have been used by Norman et al. (1980), Burkert 
& Hensler (1988), and Boss (1993) in simulations of gravi- 
tational cloud collapse and by Yanagita & Kaneko (1995) in 
their study of convection. 

We use a modified donor-cell scheme in which a planar 
function is used to approximate the gas density at each lat- 
tice site (see Appendix A). While the donor-cell scheme is 
an Eulerian technique with a fixed grid, it is most conve- 
niently defined by imagining that on each time step, the i th 
cell is transported a distance u l At, where u % is the fluid 
velocity. Each cell donates the portion of the fluid vari- 
able which overlaps its downstream neighbor. The advection 
scheme was constructed to minimize the effects of numerical 
diffusion in regions of the flow that are "ballistic" (which 
accounts for a large fraction of the life of a fluid parcel). 
In such regions numerical diffusion may be approximated 
by the d 4 /d 4 operator, while near discontinuities, it is bet- 
ter represented by the d 2 /d 2 operator, as shown in App. 
A. The effect of this prescription is that strong discontinu- 
ities quickly diffuse as Al <x At 1 ^ 2 until the discontinuity 
is spread typically over a few lattice sites, but then diffuses 
more slowly as Al <x At 1 ^ 4 once it has been softened. The 
details of the method and its associated numerical diffu- 
sion are given in Appendix A. The advantage of this ap- 
proach is that it is more computationally efficient than full 
second-order techniques, and is substantially less diffusive 
than purely first-order methods. Standard tests of hydordy- 
namic codes (e.g. the shock tube problem) cannot be used, 
since they involve pressure. However the code does accu- 
rately conserve mass and momentum, and we have checked 
the effect of resolution between 128 2 and 512 2 , as discussed 
below. However we point out that we did not carry out a con- 
ventional resolution convergence test, in which the number 
of resolution elements is varied with all other input quanti- 
ties held constant. Instead, by varying the domain size L, we 
are increasing the resolution, but are also changing charac- 
teristic scales, for example the ratio of the size corresponding 
to the peak in the energy spectrum of initial conditions to 
the size of the domain, or the total energy injection by the 
winds over the entire domain. Nevertheless, it will be seen 
below that the average values of all the statistical quantities 
of interest, including velocity probability distributions, are 
fairly stable over this factor of four change in the size of the 
system at fixed grid spacing, indicating that the results are 
not sensitive to the resolution defined in this way. Most of 
the results to be discussed are for 256 2 runs. 

2.2 Star formation 

We assume that star formation is driven by gravitational 
instabilities in expanding gas shells (e.g. Elmegreen 1994, 
Comoron & Torra 1994, Whitworth et al. 1994, McCray & 
Kafatos 1987) or in compressed layers between colliding gas 
flows (e.g. Elmegreen 1993, Vishniac 1994, Whitworth et al. 
1994). Since neither pressure nor gravity are explicitly in- 
cluded in the simulations, we rely on analytic results based 
on linear stability analysis of the hydrodynamic equations to 
determine when the gas is unstable to collapse. The adoption 
of such star formation laws based on linear stability analy- 
sis or more ad hoc schemes such as the Schmidt law must 



be used even in fully hydrodynamic simulations involving 
more physical processes (see references given below), since 
the tremendous range of scales involved in the star forma- 
tion process makes it unfeasible to follow the entire collapse. 
However our star formation criterion is different from previ- 
ous work because the criterion depends on the column den- 
sity through individual structures measured along the plane 
of the simulation. 

In general, we assume that star formation proceeds 
when the growth rate of the fastest-growing mode Ui &st 
(which depends on the geometry and overall dynamics of 
the local gas distribution) exceeds a critical growth rate ui c : 

LU[ ast > UJ C . (4) 

In analytic models of star formation which consider the sta- 
bility of expanding shells, the value of cu c usually is set to 
the age of the shell with the idea being that at earlier times 
the instability would not have had enough time to develop 
(Elmegreen 1994, Comeron & Torra 1994). The existence 
of a critical growth rate also results if collapsing gas clouds 
can be disrupted by passing shells or clouds (Franco and 
Cox 1983). Franco and Cox proposed that star formation 
proceeds only when the gas free-fall time becomes less than 
a "reagitation time" due to the passage of shells. We adopt 
this interpretation of the critical growth rate and further 
set uj c to a constant throughout the simulation. This is an 
admittedly large simplification, but we do not wish to intro- 
duce any more free parameters or spatially-dependent pro- 
cesses than necessary, so that the interpretation of the re- 
sults will be as unambiguous as possible, in keeping with our 
goal of an intermediate-level representation. 

The stability criterion used for the majority of the sim- 
ulations is based on the stability of an expanding, accreting 
sheet and is derived in Appendix B. This instability condi- 
tion depends on the velocity dispersion within the shocked 
gas layer. However, the structure of the the gas velocity field 
behind shock fronts is likely to be disordered and possibly 
turbulent since nonuniformities in the pre-shock gas distri- 
bution are known to generate vorticity (Truesdell 1952, see 
Kornriech & Scalo 2000 and references therein). Since the 
internal structure behind radiative shocks in inhomogeneous 
media is not clear, we treated the velocity dispersion of the 
shocked gas as a model parameter, constant for all shocks. 
In the future, we will explore the effects of letting the the in- 
ternal shell velocity dispersion depend on the velocity of the 
ambient gas in the shock's frame of reference as suggested 
by Elmegreen (1994). 

It is assumed that the gas is distributed in a thin sheet 
perpendicular to the simulated galactic plane and that the 
radius of curvature of the shell is large compared to its thick- 
ness. In the absence of accretion or stretching, the condition 
for gravitational collapse, and thus star formation, becomes 
a threshold condition on the shell column density: 



where fj, is the gas column density of the sheet and c is the 
shell internal velocity dispersion. 

The inclusion of accretion onto the slab and expansion 
of the slab have been studied by Elmegreen (1994) and Com- 
eron & Torra (1994) in the context of a shell expanding into 
a uniform medium. Appendix B gives a derivation similar to 
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theirs, but which leaves the accretion and expansion rates 
independent. The action of accretion and expansion do not 
affect the wavenumber or mass of the fastest growing mode, 
they just retard the growth rate making collapse more dif- 
ficult. The condition for significant collapse becomes (see 
Appendix B) 

(~~c~) > ( Wc + V ' v °+ A /m)(wc + V-v ). ( 6 ) 

where V ■ vo is the rate of stretching in the shell due to 
expansion and A is the accretion rate of ambient gas onto 
the shell. This is the condition for star formation adopted 
in the simulations. In practice the stretching and accretion 
terms were found to be relatively unimportant. A physical 
argument for the slight effect of stretching was given by 
Whitworth et al. (1994). 

Our use of the shell gravitational instability criterion 
in the simulations does not imply that the shells will frag- 
ment at sizes given by Ehlerova et al. (1997) and others 
for shells expanding into smooth ambient media. This is be- 
cause in the present simulations star formation occurs at 
multiple sites within the simulation domain, and the mean 
separation of these sites is generally small. Thus a given 
shell in the present simulations expands into an extremely 
inhomogeneous density and velocity field, and its evolution 
is dominated by interactions with the other shells of various 
ages, which have in turn been deformed by previous interac- 
tions. In addition, new shells propagate, in part, within their 
parent shells, whose densities can be large relative to the av- 
erage density. In contrast, the model shells of Ehlerova et al. 
(1997) and similar papers expand into a smooth medium and 
only become unstable when their sizes are of order 1 kpc. 

Prescriptions for SF that have been used in previous 
work on kpc-scale galaxy evolution, protogalactic collapse, 
galaxy interactions, isolated disk (and LSB disk) galaxies 
and dwarf galaxies usually involve either some Schmidt-like 
assumption about the density dependence of the SFR, (e.g. 
Chiang & Prendergast 1985, Theis et al. 1992, Mihos & 
Hernquist 1994, Anderson & Burkert 1988), or some combi- 
nation of constraints including sound crossing time exceeds 
free-fall time (a Jeans criterion) negative local velocity di- 
vergence (convergent flow), or density exceeding a constant 
threshold density (usually to imitate a condition that the 
cooling time be much smaller than the dynamical time). 
Examples are Katz, Weinberg & Hernquist (1996), Vazquez- 
Semadeni, Passot & Poquet (1996), Steinmetz & Muller 
(1995), Gerritsen (1997), Lia, Casrraro & Chiosi (1998), 
Navarro & Steinmetz (1997), Mori et al. (1997). 

In contrast to all this work, our simulations assume 
that star formation only occurs in gravitationally unstable 
sheets or filaments, for which it is essentially a critical col- 
umn density through the sheet or filament that controls local 
gravitational instability by eq. 6, for a given prescribed and 
constant gas velocity dispersion within the sheet or filament 
(since effects due to shell dilatation and accretion were found 
to be negligible) . The major cpu bottleneck in this prescrip- 
tion is that a filament-finding algorithm must be used in 
order to check all the column densities (sec. 2.3.1 below). 
Major simplifications are that the velocity dispersion in fil- 
aments is a constant, so the threshold column density is 
a constant, and that the timescale for instability is given 
by the growth timescale for the fastest-growing mode, also 



assumed constant. Because we assume that the flows are 
advection-dominated with no pressure, there is no energy 
equation or cooling condition to account for. These reduc- 
tions in the level of modeling reduce the number of free pa- 
rameters, allow a more "transparent" attempt to understand 
the physical phenomena which result (which turn out to be 
complex enough without including additional processes), yet 
contain a physically-reasonable representation of a mode of 
star formation that probably occurs in real galaxies (see the 
resent summary of references in Elmegreen 1999). 

The total mass of stars formed in an unstable region is 
assumed to be given by M sta rs = fM gas , where the star for- 
mation efficiency e is assumed constant. The gas mass M 9as 
is given by the product of the gas column density through 
the shell p s h and the vertical cross-sectional area of a simu- 
lation cell hAx, where h is the assumed scale height of the 
galaxy. We feel that this estimate of the gas mass is the most 
physically plausible given our contention that the unstable 
regions reside in gas shells. An alternative estimate would 
be to simply use the gas mass in the local star-forming site; 
however, the resulting masses would be affected by numer- 
ical diffusion, which controls the shell's thickness and thus 
average density. By using the column density, we have es- 
sentially eliminated the direct effects of numerical diffusion 
from the cluster mass estimate. 

The Jeans mass is not relevant to this calculation since 
the desired mass is the total mass of stars formed in the 
simulation cell and not the mass of an individual condensa- 
tion. For a 2-D sheet, the Jeans mass is simply Mj = fiXj, 
where is the shell surface density and Aj is the Jeans 
length. Since the number of Jeans condensations Nj in a 
2D region of size I is (l/Xj) 2 , the total mass of collapsing 
gas is MjNj = /iAj(Z/Aj) 2 = fil 2 which is simply the total 
mass contained in the region and is independent of the Jeans 
mass. Still, we refer to the newly formed stars in a given cell 
as a "cluster" with the understanding that it may be several 
unresolved clusters or a partial cluster. If the simulations 
are regarded as applying to smaller regions after suitable 
scaling, for example the interiors of molecular clouds, then 
"cluster" may be interpreted as "protostar." 

By assuming for simplicity a constant power-law IMF 
with logarithmic slope F = -1.5 and upper and lower mass 
cutoffs of 0.1 Mq and 40 Mq, we estimate the number 
of massive stars formed with M > 8 Mq to be iV* = 
M a tars/10 3 Mq. Substituting the parameter values adopted 
for the present simulations Ax = 7.8 pc (see below), h 
= 100 pc, and e = 0.1 (see section 2.5), we find that 
A* = 0.8 JVfii,2i, where Am, 21 is the filament column density 
in units of 10 21 cm -2 . This is the number of massive stars 
which form per lattice site along the filament. Thus, if a fil- 
ament becomes unstable over many lattice sites, the total 
number of massive stars produced by the filament could be 
considerably larger. 

2.3 Filament properties 

2. 3. 1 Identification 

The identification of filaments, or "shells," in the simulation 
is required to determine the shell column densities and ex- 
pansion rates used in the star formation criterion (eq. 6). 
The technique employed looks for local maxima in the den- 
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sity field along multiple lines of sight. If a given cell is a local 
maximum along at least two of four possible lines (parallel 
to the x-axis, y-axis, and diagonals) it is classified as lying 
on a filament "ridge." Requiring that a cell be a local max- 
imum in all four directions would locate cloud peaks. On 
the other hand, virtually all points where the density field 
is convex are local maxima in at least one direction (the di- 
rection tangent to the local density contour). This method 
depends entirely on the coarse-graining of the density field 
provided by the discreteness of the simulation. 




Figure 1- A typical 256 2 density field from a simulation with 
stellar forcing (top). Filaments identified by the filament- finding 
routine (bottom). Plotted are line segments of length equal to the 
lattice spacing Ax at each lattice site determined to lie on a fila- 
ment ridge. The angle of each segment indicates the estimate of 
the local filament orientation. The local filament column density 
is determined by integrating the density distribution along a line 
perpendicular to the line segments shown. 

As an example, we present the application of this tech- 
nique to an idealized filament oriented along the x-axis with 
a density gradient along the length of the filament and a 
linear falloff along the y-axis. The density profile is given by 

f(x,y) = fo+ax-b\y\. (7) 

If this function is coarse-grained by overlaying boxes of 
width Ax, the average density in a box is f(i,j) = 
J J f(x,y)dxdy, where the integration limits correspond 



to the sides of the box. If we consider a box centered at 
(x,y) = (0,0), clearly the boxes above and below it on the 
y-axis will have average densities less than the center box 
making it a local maximum along y. On the other hand, the 
average density in boxes along the x-axis increase with in- 
creasing x, so the local box is not a local maximum along 
the x-axis. Thus, the classification of this simple example as 
a filament depends on the average density in the diagonal 
boxes. It is straightforward to show that the central box will 
be a maximum along the diagonals, signaling the presence 
of a filament, if and only if \a\ < 0.756. Thus, if the density 
gradient along the filament a is too steep, the "filament" is 
not recognized. Similar arguments can be made for idealized 
filaments with arbitrary orientations and profiles. 

Figure 1 shows a section of the density field taken from 
a typical simulation and the corresponding filament ridges 
identified by this method. The local orientation of each fil- 
ament indicated by the directionality of the line segment is 
determined by the method discussed below. As the figure 
shows, the identified filament ridges reproduce the structure 
in the density field well. 

2.3.2 Shell orientation 

The local directionality of a shell is required to estimate its 
properties such as the column density and expansion rate. 
Attempts to estimate this direction by fitting linear func- 
tions to the local gas density distribution or computing the 
eigenvectors of the local gas inertia tensor were found to be 
computationally expensive and often gave erroneous results 
due to confusion with neighboring filaments. We found that 
estimating the subgrid position of the filament ridge based 
on the gas density at the neighboring lattice sites, and least- 
square-fitting a line to those neighboring positions, gave sat- 
isfactory results with low computational overhead. 

The method is as follows. For each direction in which 
the local lattice site is found to be a local maximum, the 
position of the local maximum along that line is estimated 
through parabolic interpolation. The maximum along a give 
line is given by: 

.Ax p(x" + n) - p(x - n) 

x w xq — n — r — r (8) 

2 p(x + n) + p(x - n) - 2p(x Q ) 

where n is a unit vector defining the line of interest. The 
position of the filament ridge is then defined by the average 
of these position estimates and the local filament orienta- 
tion is obtained by least-square-fitting a line to the local 
and neighboring shell position estimates. We label the di- 
rections parallel to and perpendicular to this line by unit 
vectors nil and tii for the discussion in the next two sec- 
tions. Application of this method to test density rings with 
a range of widths and radii produced errors in the estimate 
of the position angle less than 5°. 

2.3.3 Column densities 

The column density of the shell is estimated by integrat- 
ing the density distribution along a line perpendicular to 
the shell (i.e. parallel to n±). The integration is initiated 
at the shell maximum and proceeds down both sides of the 
shell until local minima in the density are encountered or a 
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maximum integration length is reached. The planar approx- 
imating functions used in the donor-cell advection scheme 
are used in the integration. Thus, the shell column density 
may be expressed as 



N(x) 



1 



p(x + sh±) ds 



(9) 



where p(x) is given by Eq. 1 and m p is the average particle 
mass. For most of the simulations, the maximum integra- 
tion distance was set to Z max = lOAx. Since the integration 
limits depend on the density structure in the simulation, 
the implementation of the algorithm which calculates the 
column densities cannot be fully vectorized and thus is rel- 
atively expensive. The requirement of a maximum distance 
in the column density integration serves to reduce computa- 
tion time. By running several simulations with varying Z max , 
we found that this limit does not affect our results. 



2.3.4 Shell expansion rate 

The time scale for the shell stretching is approximately 
i"str = [V|| • S] , where Vy represents the gradient of the 
component of the velocity field along the shell (parallel to 
nil). This gradient operator is a directional derivative and is 
defined as 
d 

V ll ' " = 'H - Pll ' v(x + Sftll) «=0 . (10) 
OS 

We estimate this derivative using a first order difference 
method 



V|| ■ V : 



(s + As) — v\\ (s — As 
2As 



(11) 



where vu = 7T.ii ■ v. Setting As = Ax and using a linear 
interpolation scheme to estimate the velocities, the following 
expression for the directional derivative results 

Vii • v » n x n y K +1 ' 3+1 - vt 1 '*- 1 ) + 



(l-n x )n y (^ +1 -v^- 1 ) + 
for n x n y > and 

Vii • v » n^nJv\-^ +1 - v l + 1 ' 3 ~ 1 ) + 



(12) 



1,7 + 1 1,7 — 1 

— Vu 



(1 — n x )n y (v 
^(l-ny)^- 1 ^ -v'* 1 '*) 



+ 



(13) 



for n x n y < where we have dropped the "||" symbol from 
the unit vector components for notational convenience. 

2.4 Stellar energy injection 

Supernovae provide the primary energy source by which 
stars drive bulk motions of the interstellar gas (Abbott 1982, 
Castor 1993). Due to the correlated positions of supernovae 
in a cluster, the net effect is a "cluster wind." Leitherer et 
al. (1993) showed that the net power and momentum of the 
combined winds and supernovae ejecta of a star cluster re- 
main roughly constant over the first 10 7 years and declines 
afterward (see, however, Oey & Smedley 1998). They found 
this result to be largely independent of the the model pa- 
rameters including the metallicity and the detailed form of 
the IMF. 



We model the energy injection from a star cluster as a 
cluster wind (conceptually viewed as the product of multi- 
ple supernova events) with constant momentum and power 
which lasts r w years. The total energy released over the life- 
time of the cluster is simply E to t = N,Esn, where N* is 
the number of massive stars in the cluster which become su- 
pernovae and Esn is the energy per supernova. We assume 
that the transfer of the cluster wind energy into mechanical 
kinetic energy of the gas within the local star-forming lat- 
tice site has an efficiency / mcc h. This efficiency parameter 
is introduced to parameterize the effects of the unresolved 
physics related to heating and radiative losses during the 
formation of the expanding shell. Thus, the kinetic energy 
of the local gas is given by E w = / mec h-Etot and its average 
velocity is v w = (2_B m /m cc n) 1//2 , where m ce il is the mass in 
the local star forming cell. Substituting the above expression 
for Etot, the expression for the gas velocity becomes 



2/mcchA?».EsN \ 1//2 



Tlccll 



(14) 



By substituting the expression for the number of mas- 
sive stars formed in a cluster (see section 2.2.1), N„ = 
em ce ii/lO 3 M0, we find that typical values of this velocity 
are given by 

j) (4^Wn,si1 1/2 



= 40 km s 



0. 



0.1 



(15) 



Thus 40 km sec -1 would be the velocity of the gas at a 
typical star forming forming lattice site if it received a frac- 
tion / mcc h = 0.1 of the total energy produced by the su- 
pernovae in a typical star cluster. Our standard value of 
fmechEsN = 10 50 erg is close to the value found in de- 
tailed 1-dimensional simulations of supernova remnants by 
Thornton et al. (1998). However our wind velocities are 
much smaller than usually estimated (~1000 km s~ ) be- 
cause, in our simulations, the momentum input has to ac- 
celerate a mass of gas in a cylinder whose height is equal 
to the adopted scale height. This small wind velocity is, on 
the other hand, computationally expedient because larger 
wind velocities would require smaller time steps. Besides, 
our wind velocity is only applied at distances from the clus- 
ter equal to half a grid size (3.9 pc) and usually occurs within 
high-density regions, so the expected wind speeds should not 
be as large as 1000 km s -1 . 

Computationally, the cluster wind is modeled by redis- 
tributing the mass in the local star forming site into the 
site's eight nearest neighbors over a time interval r w . The 
redistributed mass is given velocity u w directed radially away 
from the cluster. After the initial mass and momentum de- 
position, the dynamics are governed by the donor-cell ad- 
vection prescription. 

2.5 Shell expansion 

As the cluster wind encounters the surrounding gas a com- 
pressed layer forms and expands away from the star forma- 
tion site, conserving momentum. For a uniform ambient gas 
distribution, shells driven by internal winds are expected to 
follow the R s h <x t 2,/3 law predicted by extending Avedis- 
ova's momentum-conserving solution (Avedisova 1972) to 2- 
D. After the cluster winds subside, the shell would evolve as 
Rsh oc t 1 ^ 3 as in the 2-D zero-pressure snowplow solution. 
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Neither of these expansion laws are expected to hold 
for true wind-blown shells. In each case, the two dimen- 
sionality of the simulations changes the expansion power- 
law exponents. The zero-pressure snowplow solution in 3-D 
is -Rah oc i 1//4 , and Avedisova's 3-D momentum conserving 
solution has an expansion law of t 1 ' 2 . Also, Weaver et al. 
(1977) demonstrated that the hot interior of the shells can 
settle into pressure equilibrium and that this pressure leads 
to a faster expansion oc t 3 ^ 5 than Avedisova's solution. 
The reason that the shells in the present simulations are ex- 
pected to follow Avedisova's solution and not Weaver et al. 
is that we do not model pressure effects. However, we do 
not feel that the precise form of the expansion law will have 
a qualitative effect on the global statistical properties and 
organizational behavior in which we are interested. Thus, 
the added computational expense of solving an energy equa- 
tion to better reproduce the shell expansion law does not 
seem warranted. Also, we note that the t 2//3 scaling expected 
in the present simulations is intermediate between the 3-D 
Weaver solution (t 6 ) and the 2-D analog (R s h oc t 3 ^ 4 ) which 
we derived using simple dimensional analysis arguments. 

These scaling solutions only apply to shells expanding 
into a uniform medium, but in the present simulations the 
gas becomes very inhomogeneous, and, besides, the shells 
usually begin their expansion within a "parent shell" and 
so even the local ambient density field is anisotropic. For 
these reasons, none of these idealized scaling solutions are 
realized, except perhaps in a statistical sense at early times. 
Rather, the shells are often stretched and distorted as they 
interact with the surrounding gas density and velocity fields. 

2.6 Normalizing values 

Since we are primarily interested in modeling the dynam- 
ics of the ISM in disk galaxies, we choose fiducial time and 
length scales appropriate to that problem. However, we note 
that through renormalization, the model could also be ap- 
plied to smaller-scale, lower-momentum interactions driven 
by low-mass protostellar winds within a molecular cloud, 
such as the model suggested by Norman & Silk (1980), or 
to megaparsec-sized scales in which galactic winds provide 
the momentum forcing, such as the propagating galaxy for- 
mation model of Ostriker and Cowie (1981). 

We set the lattice spacing to Aa; = 7.8 pc giving sim- 
ulation sizes of L = 1 kpc, L — 2 kpc, and L — 4 kpc for 
the 128 2 , 256 2 , and 512 2 grids, respectively. This choice is 
somewhat arbitrary but represents a compromise between 
choosing a large enough domain to contain many superbub- 
bles (which are typically observed to have sizes from 50 to 
100 pc but vary from less that 10 pc to more than a kilopar- 
sec (Heiles 1979, Chu & Kennicutt 1994, Yang et al. 1995)) 
and one small enough such that the absence of galactic dif- 
ferential rotation in the simulation is not too severe and 
such that the shells are resolved by many lattice sites. In the 
solar neighborhood, the shearing time of a kiloparsec-sized 
region is roughly r s h car = L(dV/dR)~ 1 ~ 2 x 10 8 years (Mi- 
halas & Binney 1981). In comparison, the kiloparsec cross- 
ing time for the fluid motions in the present simulation is 
to = 2.5 x 10 7 years (see below). While the effects of differen- 
tial rotation can have important consequences for large-scale 
fluid dynamics (e.g. Palous et al. 1994 and Elmegreen 1994 
study its effects on expanding supershells) , we wish to inves- 



tigate the effects of advection on the shells, decoupled from 
other processes. 

We adopt an average gas density of no = 1 cm -3 cor- 
responding to typical values in galaxies on kiloparsec scales 
and a disk scale height h = 100 pc, although one series of 
models was run with a variety of values of no. With this 
choice, the mass in a cell with the average density no is 
nicoii = 190 Mq. The average gas column density AT ce ii in 
the plane of the simulation through a cell with density no is 
iV ce ii = n Ax = 2.4 x 10 19 cm~ 2 n (Aa;/7.8pc). 

The cluster wind speed provides a natural character- 
istic velocity for the simulations. In the previous section a 
computationally appropriate wind velocity was found to be 
v w — 40kms~ 1 . We define a characteristic time scale as the 
kiloparsec crossing time of the wind to = lkpc/uo = 2.5 x 10 7 
years. 

We note that the two-dimensionality of the simulation 
poses several problems in assigning physical units. First, real 
wind-blown shells are initially much smaller than the galaxy 
scale height and expand into a three-dimensional medium. 
Since the mass in each simulation cell is assumed to corre- 
spond to a column of gas with a height equal to the galaxy 
scale height h and mass Ax 2 hpo, then the velocities of the 
expanding shells would be vastly underestimated initially if 
they are imagined to sweep up all the gas in the local cell 
(column). However, once the shell becomes comparable in 
size to the galaxy scale height, the simulation velocities will 
be in reasonable agreement with true shell velocities since 
the affected mass is R 2 h hpo ~ h 3 po- In addition, as discussed 
above, the expected scaling of the shell expansion laws de- 
pends on the spatial dimension. 



2.7 Reynolds number 

The Reynolds number Re measures the relative importance 
of the inertial (advective) to viscous or diffusive terms and 
may be estimated by comparing the time scales associated 
with each process. The inertial time is simply Ti ncrt i a i = L/v 
where v and L are characteristic velocities and sizes of a re- 
gion of the flow. The only diffusive terms present arise from 
numerical diffusion of the donor-cell scheme. Near discon- 
tinuities, the numerical diffusion has the form vid^ where 
v% — 0.5(1 — u/u s )u/u & , where u s = Ax/ At is the maxi- 
mum stable grid velocity. The corresponding time scale is 
tv,2 ~ L 2 /v2 = 2L 2 /(1 — u/u s )u/u s . Thus, the Reynolds 
number is Re ~ TV^/Tjnertiai ~ 2L, where we have jus- 
tifiably assumed that typical grid velocities are, on aver- 
age, considerably less than unity (in units of v w ). Thus, for 
128 2 and 256 2 simulations, the Reynolds number is about 
Re — 256 and 512, respectively, indicating that diffusion is 
relatively unimportant on the largest scales of the simula- 
tions. By equating the diffusive and inertial times, we find 
that diffusion becomes important only on the smallest scales 
Id ~ Aa;. The characteristic thicknesses of filaments in our 
simulations are a few times 

Away from discontinuities, the numerical diffusion takes 
the form -v A d 4 /dx 4 with v A = [3(1 - (u/u g ) 3 ) - 6(1 - 
u/u s )u/u s }u/u s . In this case, the "super" diffusion time is 
approximately t v ,4 ~ I 4 /va ~ I 4 /(3u/u g ) for small veloc- 
ities. The resulting Reynolds number Re fa l 3 /3 is Re w 
7 x 10 5 and 6 x 10 6 for 128 2 and 256 2 lattices respectively. 
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2.8 Initial conditions 

The velocity field is initialized as a Gaussian random field, 
meaning that the spatial Fourier modes have random uncor- 
rected phases with Gaussian-distributed amplitudes depen- 
dent only on the norm of the Fourier wavenumber k = |fe|. 
These conditions guarantee that the probability distribution 
function of a given mode is initially Gaussian. Two forms of 
the assumed initial energy spectrum are used: a scale-free 
spectrum of the form E(k) oc k~ a and an energy spectrum 
which peaks at a preferred scale ko, E(k) oc k 4 exp( — k 2 /kg). 
The standard deviation of the initial velocity field is set near 
the characteristic velocity vq to minimize the transient time. 
The value of the initial variance is found to not affect the 
post-transient behavior. No stars are initially present. For 
the majority of the runs, the initial density field is assumed 
to be uniform. Fluid advection, however, quickly causes den- 
sity inhomogeneities to develop which reflect the form of the 
initial velocity fluctuations. This is similar to the manner 
in which Elmegreen et al. (1995) initialized simulations of 
shock-turbulence interactions. 

2.9 Boundary Conditions 

Periodic boundary conditions were used for all the simula- 
tions presented in this work, so the simulations take place 
on a torus. Using this choice, the simulation domain is in- 
tended to represent a portion of a more extended region of 
space with statistically similar dynamics. While this choice 
is widely used in spatial simulations, it can introduce artifi- 
cial correlations on scales comparable to the domain size. 

In the future, we plan on implementing the code on a 
differentially rotating coordinate system in which concentric 
rings of lattice sites are allowed to slip past each other to 
mimic the effects of galactic differential rotation (e.g. Seiden 
& Gerola 1982, Palous et al. 1994). However at present, our 
main interest is in isolating the behavior of wind-advection 
interactions. 



3 RESULTS: SPATIAL DISTRIBUTION 
3.1 Model parameters 

Table 1 lists the parameter values and a few global statis- 
tics for the simulations with star formation. The models are 
organized such that each series represents the variation of 
a parameter relative to a standard model (indicated by an 
asterisk in the table). For runs A, the domain size was varied 
from 64 2 to 512 2 corresponding to linear scales from 0.5 kpc 
to 4 kpc. We define the standard model as the 256 2 simu- 
lation. Notice that all the statistical quantities listed show 
little variation for domain sizes larger than 128 2 . As pointed 
out in section 2.1, the variation of domain size is not equiva- 
lent to conventional resolution convergence tests because the 
resolution is not being varied independently of other phys- 
ical parameters. However, the stability of the quantities in 
Table 1 does show that the results are not affected by reso- 
lution in the sense of varying the domain size at fixed grid 
spacing. 

In runs C, the parameter varied is the internal gas veloc- 
ity dispersion of shells c s h which determines the shell star 
formation column density threshold iV s f. In the standard 



model this parameter value is set to c s h = 1 km s _1 which 
corresponds to iVsf,2l = 1.0. The models labelled with a 
dagger refer to runs that use the same parameters as the 
standard model except that the domain size is 128 2 rather 
than 256 2 . 

In runs D, the time delay Td between the initiation of 
local gravitational instability and the onset of the cluster 
winds was varied, with the value for the standard model 
being 2.5 x 10 6 years. As pointed out by the referee, this 
adopted standard value is probbly too small. The shell 
growth is dominated by the lifetimes of the most massive 
stars (Leitherer et al. 1993, Oey 1996), which is a factor 
of 1.5 - 2 times larger than the standard time delay value 
adopted here. If cluster winds can be driven by the collec- 
tive action of protostellar winds, instead of multiple SNe, 
the onset of the cluster wind should roughly coincide with 
the formation of the first stars. However, Td also includes the 
gravitational collapse time, which depends on the density of 
the gas from which the cluster forms. From these consider- 
ations, we think it is likely that the models listed in Table 1 
with Td — 0.5 — 2 x 10 7 yr may be the most realistic. These 
cases are discussed in sec. 3.2.4 and 4.1.3 below. 

In runs E, the average energy released per supernovae 
-Esn was varied around the standard value of 10 51 ergs. Runs 
F vary the average density no. Within some of the series, the 
effects of varying the simulation size and the power-law slope 
of the initial velocity power spectrum were investigated. We 
also studied a series of "decay" runs in which star formation 
was "turned off" to study the organization and kinematics 
of the gas in the absence of stellar forcing. 

Each of the simulations were integrated for 2 Gyrs. No- 
tice that the integration time (about 40 crossing times for 
the 2 kpc simulations) is much larger than in most previously 
published kpc-scale simulations, because we were especially 
concerned with the possibility of long transient times and 
wanted to search for any relatively large time-scale correla- 
tions induced by spatial self-organization. Such long integra- 
tion times were possible because of the intermediate level of 
representation of our model (e.g. self-gravity only included 
as a local instability prescription) and the use of a first-order 
differencing scheme. 

The physical motivations for our choice of normalizing 
parameters were given in section 2.5. Here, we briefly list 
these values. We set the star formation efficiency to e = 0.1 
and assume that one massive star forms per total 1000 Mq 
of stars formed. An initial average density of no = 1 cm -1 
is chosen for all models except for series F. Assuming that 
the fraction of the total energy released by supernova which 
is transferred to kinetic energy of the shell is / moc h = 0.1 
(consistent with Thornton et al. 1998), the characteristic 
velocity and kiloparsec crossing time of the wind are vo = 40 
km s _1 and to = 2.5 x 10 7 years respectively. The duration 
of a typical simulation is 2 Gyr or 40 grid crossing times 
for the 256 2 simulations. Assuming a galactic scale height 
of 100 pc, the mass in a cell with this average density is 
m s f = 1.7 x 10 4 Mq. A flat energy spectrum of the initial 
velocity field with equal power per logarithmic wavenumber 
bin was chosen for most of the runs, but some Ey, oc fc~ 2 
initial power spectra were also studied. 

Gas consumption by stars is not included in the models 
presented. Gas consumption was studied, but was found to 
simply lead to a gradual decline in the global star formation 
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Table 1. The wind duration is t w = 10 Myr for all runs. The columns are (2) linear grid size, (3) time delay in Myrs, (4) average 
gas density, (5) Energy per supernovae in 10 51 ergs, (6) shell internal velocity dispersion in km s —1 , (7) logarithmic slope of the 
energy spectrum of the initial velocity field. Statistics in columns 8-18 are based on the last half of the simulations: (8) standard 
deviation of the density field, (9) average number of density peaks along a line of sight, (10) mass-weighted gas velocity dispersion, 
(11) velocity dispersion of the shells, (12) average gas kinetic energy, (13) average number of massive stars per kpc~ 2 , (14) average 
number of star clusters per kpc" 2 , (15) average energy injection rate, (16) relative width of the SFR histogram, where b = N t /N*. 
The (T) and (*) symbols indicate the standard model at resolutions of 128 2 and 256 2 respectively. 
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rate and gas velocity dispersion, as expected. Since this de- 
cline occurs only over large secular time scales, we turn off 
stellar gas consumption for the runs presented here in order 
to produce simulations whose temporal behavior is nearly 
statistically stationary, so that time averages can be rea- 
sonably defined, and so that simulation data samples from 
different times can be combined to decrease noise in sta- 
tistical estimates for quantities like the velocity and mass 
probability distributions. However it is important to keep in 
mind that results that are found to depend on the average 
SFR for different models apply to a single model at different 
times if gas consumption were included, since in that case 
the average SFR decreases with time. Models including gas 
consumption are described in Chappell (1997). 

Since the filaments in the 2-D simulations are meant to 
represent expanding and interacting shells in 3-D, we will 
use the terms "filaments" and "shells" interchangeably in 
referring to these density structures. However, it should be 
remembered that in 3-D the morphology will likely be a 
mixture of shell-like (i.e. two-dimensional) and filament-like 
(i.e. one-dimensional) structures. 



3.2 Spatial distributions of the stars and gas 

3.2.1 Standard model 

Figure 2 shows the logarithm of the gas surface density (grey 
scale image in each panel) and recent star formation activity 
(right panels) at two times separated by 7 x 10 7 yrs for the 
standard model. The filled circles, open circles, and crosses 
indicate the birth positions of clusters with ages between 
zero and 10 7 years, 10 7 and 5 x 10 7 years, and 5 x 10 7 and 10 8 
years respectively. The clusters in the youngest age bin have 
active winds while the winds from the oldest clusters shown 
have long since subsided. The symbols mark the positions 
of the clusters at the time when their winds were active 
in order to show how the spatial distribution of the star 
formation sites evolve in time. The change in the cluster 
positions (the clusters inherit the velocity of their parent 
gas clouds) between the time of their birth and the time 
that the density field is shown is typically less than a lattice 
spacing for the youngest clusters, but can range up to a few 
hundred parsecs for the oldest (10 s year) clusters. 

The density fields shown in Fig. 2 exhibit a network or 
web of shells or filaments which are distributed nonuniformly 
over the grid. Large "voids" and concentrations of shells can 
be seen at both times. The gas distribution qualitatively re- 
sembles the aperture synthesis H I map of the LMC pre- 
sented by Kim et al. (1998). The most recent star forma- 
tion events (filled circles) generally trace the large-scale gas 
distribution. However, on small scales, winds from all but 
the very youngest clusters have excavated low-density cav- 
ities around the clusters and initiated expanding gas shells 
(see, for example, the clusters and surrounding shells in the 
boxed region at time t = 1.99 Gyrs.) We find that the fate 
of these shells depends largely on the structure of the sur- 
rounding gas density and velocity fields. In some instances 
the shell merges with the ambient gas but never reaches the 
star formation threshold, and thus does not directly trigger 
the formation of more clusters. In other cases, shell mergers 
lead to the formation of lone clusters as evidenced by the 
presence of isolated clusters in the right hand panels. Much 



more common, however, is the tendency for star clusters to 
form in large groupings. 

Figure 3 shows the time evolution of the boxed region 
of the density field between the two times shown in Fig. 2. 
In the lower region of the the initial panel, two young star 
clusters have just initiated two expanding shells. These shells 
collide at At = lOMyrs, triggering the formation of two new 
clusters. By At = 30 Myrs, the shells have merged to form a 
single "supershell" which is expanding toward the filament 
in the center of the panel. At time At = 40 Myrs, the merg- 
ing of these two density structures triggers a new round of 
star formation in the center of the panel and in the lower left 
corner. The expanding shells from these two star formation 
events leads to a cascade of multiple propagating star forma- 
tion events, until, at time At = 70 Myrs, a series of nested 
shells is apparent. The organization of propagation events in 
this model is controlled to a large degree by the surrounding 
gas distribution, which is determined by older propagation 
events. Thus, propagating star formation in this model does 
not evolve as a uniform progression through a cloud as orig- 
inally suggested in the "burning cigar" model of Elmegreen 
& Lada (1977), nor does it lend itself to descriptions in- 
volving the instabilities of expanding shells in a uniform 
medium (Comeron & Torra 1994, Elmegreen 1994). Rather, 
even on the relatively small scales of individual propagating 
events, the process appears as a complex cascade mediated 
by multiple interacting shells. The specific ionization-shock 
fron model of Elmegreen & Lada (1977) is of course not re- 
futed by the present calculations, since their model involves 
scales which we cannot resolve; we only contrast our result 
with their model in a generic sense. 

The emergence of large-scale organization in the star 
formation activity is apparent in Fig. 2. The large coher- 
ent groups of stellar clusters is unexpected. These groups 
stimulate new star formation in neighboring regions, causing 
the locus of star formation activity to "wander" throughout 
the simulation. However, this apparent motion is complex in 
the sense that it depends on the organization of the ambi- 
ent gas, which evolves through the action of fluid advection 
and the collective action of many stellar generations. Thus, 
the organization and progression of the large-scale star for- 
mation activity appear not to occur as a single propagat- 
ing star formation wave in which stars form consecutively 
behind an advancing front, but rather as the organization 
of many smaller-scale propagation events which interact in 
non-trivial ways. 

For example, the winds from the old group of clusters 
indicated by the x's in the lower-right corner of Fig. 2 at 
time t = 1.92 Gyrs has created a large, low density cavity. 
The more recent star formation activity indicated by the 
closed and open circles to the lower left of the boxed 
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Figure 2. Logarithm of the gas density (left panels) and recent star formation activity (right panels) at two times separated 
by 7 x f 7 yr for the standard model. Filled circles, open circles, crosses: birth positions of clusters with ages between zero 
and 10 7 yr, 10 7 and 5 x 10 7 yr, and 5 x 10 7 and 10 8 yr, respectively. 
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Figure 3. Time evolution of the density field for the boxed region indicated in Fig. 2. The sequence begins at time t=1.92 Gyr 
(top panel of Fig. 2) and progresses in increments of 10 Myr. The final panel at t=1.99 Gyr corresponds to the lower panels 
of Fig. 2. 
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region lies at the left edge of this cavity. The "triggering" 
of this more recent star formation region was a result of the 
shells driven by the older star clusters at the the center of 
the this cavity as well as the clusters to the upper left and 
at the top of the grid. 

The two-point correlation functions of stars in the sim- 
ulations are discussed elsewhere (Scalo & Chappell 1998a). 
Their form is generally power-law, but the logarithmic slope 
depends on the age of the stars selected and on the average 
level of star formation in the models. 



3.2.2 Decay runs 

The filamentary structure of the gas in these simulations 
is a result of advection in the presence of high compress- 
ibility. Even in the absence of forcing by cluster winds, ini- 
tial fluctuations in the gas velocity field lead to filamentary 
structure. Figure 4 shows the evolution of the density fields 
for two 256 2 simulations in the absence of star formation. 
The initial energy spectra of the left and right panels are of 
the form Ey_ oc k 4 exp( — k 2 /k(,) peaking at ko = 8 L _1 and 
oc k, respectively. At early times, soon after shock for- 
mation, the density field reflects the structure of the initial 
velocity field, exhibiting filamentary structures of advected 
gas at the shock fronts with the spatial wavenumber of the 
filaments being approximately equal to the initial wavenum- 
ber ko in the left panels. As the shocks collide, they conserve 
mass and momentum and merge due to the high compress- 
ibility of the gas. Thus, the number of shocks (and equiv- 
alently filaments) decreases with time and the correlation 
length of the density and velocity fields increase. The inter- 
section of oblique filaments creates density clumps. These 
clumps are similar in nature to the clumps that form at the 
intersections of filaments in large-scale cosmological simula- 
tions. At late times, only a few dense clumps and filaments 
remain. While the filamentary structure of the gas is a direct 
consequence of the high compressibility of the simulations, 
comparison of Fig. 2 and Fig. 4 clearly shows that star for- 
mation feedback influences the organization of the filaments. 

We find that the kinetic energy decays as t~ x for all 
the decay runs that have an initial energy spectrum peaked 
at a length scale fc" 1 <C L, so that many shock interac- 
tions occur. This decay looks similar to the results found 
in more detailed 3D fully hydrodynamic and MHD simula- 
tions (MacLow et al. 1998, Stone et al. 1998). This simi- 
larity suggests that the insensitivity of the 3D decay to the 
strength of the magnetic field or the Mach number can be 
understood as the result of the dominance of fluid advec- 
tion in controlling the dynamics. A detailed discussion of 
the kinetic energy decay and power spectrum of the decay 
simulations and comparison with various analytic results for 
Burgers turbulence is presented elsewhere, since our major 
theme here is the influence of star formation. The proba- 
bility distribution (pdf) of the density field for the decay 
simulations was presented in Scalo et al. (1998), where it 
was shown that this density pdf is a power law at large den- 
sities because of the effective polytropic exponent being less 
than unity (zero here), in contrast to isothermal simulations, 
which are in a sense a singular case, giving lognormal den- 
sity pdfs. The reader is referred to Scalo et al. (1998) and 
Passot & Vazquez-Semadeni (1999) for details. 



3.2.3 Star-formation threshold 

The gas column density threshold for star formation in this 
series of models (series C in Table 1) is found to control both 
the global star formation rate and the spatial organization 
of the stars and filaments. Figure 5 shows the organization 
of the gas filaments (left) and star formation sites (right) for 
three values of the star formation threshold. As the thresh- 
old inceases, the number of filaments and the spatial density 
of star clusters decreases and the characteristic size of the 
voids between filament increases. In the upper panels which 
show the largest values of the star formation threshold, only 
a few young clusters with active winds are prsent at any 
given time. The shells expand to large radii before interact- 
ing and the characteristic time for interactions is correspond- 
ingly long. Much of the shell's life is spent in a state of free 
expansion in which the driving winds have long since sub- 
sided. At the lowest star formation threshold, young shells 
rapidly interact, triggering vigorous star formation activity 
and new expanding shells. The characteristic scale of recog- 
nizable shells is much smaller and they often interact while 
still being driven by stellar winds. 

The shell size is largely controlled by mass conserva- 
tion. The threshold parameter sets the upper limit on the 
gas column density of shells; above that vaue, star forma- 
tion is triggered and the shell is disrupted by winds from the 
young embedded stars. For a shell expanding into a uniform 
medium, the higher threshold column density is reached at 
a larger shell radius and at a later time than for a lower 
threshold value. In the present simulations in which most of 
the mass resides in shells, mass conservation requires that 
if each filament is on the average more massive, then their 
number density must be lower. Thus, as the star formation 
threshold is increased, the characteristic time between suc- 
cessive propagation events decreases and the shells are on 
average larger and older. 

At the largest star formation rates investigated in this 
series, the gas is highly agitated and the filaments appear 
shredded. Thus, a typical fluid parcel is not able to decel- 
erate significantly before being swept-up in another shell. 
Fluid parcels are frequently driven by more than one wind 
source in these simulations. 

3.2.4 Time delay 

As discussed in section 2, once the gas in a simulation cell 
satisfies the star formation criterion, a fraction (90% in the 
present simulations) of the gas in that cell is removed from 
the grid to form a bound cloud which no longer interacts 
with the surrounding gas flow. This cloud inherits the space 
velocity of the gas in the local cell and moves as a colli- 
sionless particle. The time delay parameter Td studied in 
this section is the time between the onset of the collapse of 
the gas (and formation of the bound cloud) and the onset 
of stellar winds. Over the lifetime of the cluster, the mass 
of the bound cloud is returned to the ISM as the cluster 
wind. Thus, larger time delays Td allow the gas around newly 
formed proto-clusters to evolve longer before the onset of the 
cluster winds. 

The density distribution and star formation sites for a 
model with time delay Td = 10 Myr is shown in Fig. 6 for 
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Figure 4. Evolution of the density field for two 256 2 simulations in the absence of star formation ("decay runs"). Left panels: 
Initial energy spectrum Ek ~ k exp(—k 2 /k„) with k = 8I/ _1 . Right panels: Scale free initial energy spectrum Eu ~ fc. 
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Figure 5. Typical snapshots of gas density (left) and birthsite (right) distributions for simulations in which the column density 
threshold for star formation is decreased (from top to bottom). N s f ; 2i is the critical filament column density at which star 
formation was assumed to occur, in units of 10 21 cm -2 . 
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Figure 6. Gas density distribution (left) and star formation sites (right) at three times for a model in which the time delay is 
10 Myr, compared to 2.5 Myr for the standard model (cf. Fig. 2 and middle panels of Fig. 5). Overall star formation rate and 
clustering level are enhanced. 
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three times. As the time delay has been increased relative 
to the standard model, the cluster formation rate is larger. 
Accompanying the increased SFR is a higher presence of 
"chains" of star clusters (compare Fig. 6 where Td = 10 Myrs 
with the standard model in Figure 2 where Td = 2.5 Myrs). 
These stellar chains trace the underlying gas distribution of 
star-forming filaments. 

Since the winds from the young stars in a filament dis- 
rupt the gas distribution in the filament, a larger time delay 
increases the likelihood that the column density of a fila- 
ment, which is sweeping up ambient gas, will reach the star 
formation threshold before it is disrupted. Thus, more star 
clusters will form on average in a given filament. For simi- 
lar reasons, the longer time delays increase the chance that, 
when the winds from a young cluster compress the gas along 
the filament, the star formation column density threshold 
will be reached sooner leading to more star clusters in a 
given filament. In other words, longer time delays allow fila- 
ments to accrete more mass (before being disrupted by the 
winds from the first stellar cluster to form), which increases 
the probability that more stars clusters will form "sponta- 
neously" along the filament and which decreases the length 
scale along the filament for "triggered" star formation as 
the stellar winds compress the gas along the filament. This 
scenario also explains the higher global SFR since each star- 
forming filament produces more star clusters on average. 

Clearly, the regularity of these star cluster chains de- 
rives from the uniform gas distribution along the filaments. 
The inclusion of subgrid processes would undoubtedly gen- 
erate finer-scale structure within the filaments which would 
lead to less regular stellar distributions. However, our quali- 
tative finding that longer delays between star formation and 
stellar energy injection produce a higher global SFR with 
larger spatio-temporal correlations should remain valid. 

If the cluster winds are driven by the collective action 
of the stellar winds, then small time delays may be appropi- 
rate, since stellar winds set in during the protostellar phase. 
However, if the cluster winds are driven by collective SNe, 
as in the standard picture, the longer time delays are appro- 
priate, since the evolulutionary time required to produce a 
SN varies from about 3 x 10 6 yr for the most massive stars 
to about 3 x 10 7 yr for the least massive stars believed to be 
capable of SN explosions. Overall, we believe that the long 
time delay simulations are most likely to represent the real 
situation, although only a few such simulations are presented 
here. 

3.2.5 Discussion 

The filamentary structure, observed in this simulation as a 
result of high compressibility and forcing by stellar winds, is 
found in many other simulations of the ISM in which stars 
drive gas motions. The simulations of Chiang & Prender- 
gast (1985) and Chiang & Bregman (1988), which treat the 
stars as a continuous fluid that are born at a rate propor- 
tional to some power of the local gas density and include pa- 
rameterized power-law cooling and stellar heating, develop 
a network of filaments near pressure equilibrium. Due to 
the regularity of these structures, they predict that in three 
dimensions the voids would be roughly spherical in shape 
and that the filaments seen in their 2-D simulation would 
correspond to shells or sheets. More recent work by Rosen 



et al. (1993) extended their approach by using an integra- 
tion scheme with less numerical diffusion and found that the 
resulting filamentary network is less symmetrical and that 
filaments with a range of densities and temperatures emerge. 

The star formation-powered turbulence simulations of 
Vazquez- Semadeni et al. (1994), Passot et al. (1996) and 
Vazquez- Semadeni et al. (1997), which include more physics 
(self-gravity, magnetic fields, expicit heating and cooling, ro- 
tation) also exhibit filamentary structures, although they are 
transient and less regular than the structures found by Chi- 
ang & Prendergast (1985). These structures tend to form at 
the interface between colliding gas streams and lead to star 
formation, similar to the picture suggested by Elmegreen 
(1993). Similar structure is found in the isothermal, non-self- 
gravitating, randomly forced MHD simulations of Padoan et 
al. (1998). Given the ubiquity of these structures in simula- 
tions that involve very different formulations of the physics 
(i.e. cooling, pressure, star formation "laws", self-gravity, 
etc.), the present simulations, which only incorporate highly 
compressible advection and driving by star formation, sug- 
gest that these latter physical effects are responsible for most 
of the filamentary structure seen in these simulations. 

Finally, we consider what the gas distribution would 
look like were it observed with non-perfect resolution. The 
gray-scale image in Fig. 7 shows the density field convolved 
with a Gaussian with FWHM of 250 pc. This linear resolu- 
tion corresponds to angular resolutions of 30', 26", and 5" 
at distances of 50 kpc, 2 Mpc, and 10 Mpc respectively. No- 
tice also that because the structure is severely underresolved 
at the resolution of the simulated observations, the gas den- 
sity actually appears to be distributed in roundish "clouds" ; 
comparison to panel a suggests how misleading conceptual 
models based on the under-resolved observations might be. 



4 EVOLUTION OF GLOBAL STATISTICS 
4.1 Global star formation rate 

4- 1.1 Standard model and effects of domain size 

Figure 8 shows the evolution of the global star formation 
rate (left-hand axis) and gas velocity dispersion (right-hand 
axis) for the standard model with the domain size ranging 
from 250 pc to 2 kpc. The star formation rate is defined 
as the number of massive stars with active winds per kpc 2 . 
The time-averaged value of the number of clusters with ac- 
tive winds per unit area is around 7 kpc -2 and is largely 
independent of the simulation size except for the smallest 
grid. Since the wind duration is t w = 10 Myr, this cluster 
surface density corresponds to a cluster formation rate of 
7 x 10" 7 kpc" 2 yr" 1 which is about a factor of three greater 
than the cluster formation rate for the local Milky Way re- 
ported by Elmegreen & Clemens (1985) of 2.5 x 10" 7 kpc" 2 
yr" 1 . We consider this to be in acceptable agreement with 
observations, considering that the aim of this work is not to 
reproduce the detailed properties of any single galaxy, but 
to investigate the dynamical and organizational aspects of 
galaxies. 

The velocity dispersion shown in the same figure is de- 
fined as the rms velocity component weighted by the gas den- 
sity a 2 p(v x ) = jjJ2i( ViWi -Wm) 2 , where wt = pif(jf'Epi) 
and N is the number of lattice sites in the grid. Thus, it is 
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Figure 7. Smoothed gas density and stellar fields for a typical time in a series C simulation. Grey scale images represent 
gas density smeared with a Gaussian filter and are identical in the two panels. White corresponds to highest densities. Top 
panel: solid contours are smoothed distribution of clusters with active winds (t < 10 7 yr) and dashed white lines show 
distribution for clusters with ages 5 x 10 7 < t < 10 8 yr. Contours in lower panel are distribution of intermediate-aged stars 
with 10 7 < t < 5 x 10 7 yr. 
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Figure 8. Global star formation rate (number of massive stars per square kpc) and gas velocity dispersion (solid and dashed 
lines, respectively) as a function of time for models in series A. The simulation sizes are (a) 512 2 , (b) 256 2 , (c) 128 2 , (d) 64 2 . 
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Figure 9. Frequency distribution of the ratio of Oh? present SFR to the mean SFR. For the last Gyr of the simulation, 
denoted by b, for the four simulations with different spatial extent corresponding to Fig. 8. The burstier behavior with 
decreasing size is manifested as a broadening of the b-distribution. 
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meant to approximate the width of an optically thin spectral 
line if the simulation were viewed by an observer in the grid 
plane with a beam which contains the entire lattice. The 
average gas velocity dispersion weighted by the density field 
is around a p — 0.12«o which translates to 5 km s" 1 . The 
velocity dispersion of the filaments for the runs in series A 
correspond to 2.2 km s _1 . Stark & Brand (1989) find that 
the cloud velocity dispersion in our galaxy is about 7 km 
s _1 for cloud masses up to 3 x 10 5 Mq and declines to 3 km 
s" 1 for larger clouds. Again, we consider these results to be 
in rough agreement. 

Fluctuations in the global statistics increase in ampli- 
tude as the system size is decreased. For L less than a kilo- 
parsec, the star formation rate becomes increasingly bursty 
with long temporal correlations present, while for L greater 
than a kiloparsec, the star formation rate and gas veloc- 
ity dispersion settle into a near equilibrium state with no 
secular variations in time. This transition from continuous 
star formation to bursty behavior is reminiscent of the tran- 
sition proposed by Gerola et al. (1980) in their stochas- 
tic self-propagating star formation (SSPSF, see Seiden & 
Gerola 1982 for a review) theory of bursting dwarf galaxies. 
In both cases, large simulations, in which many regions of 
propagating star formation are present, exhibit near equi- 
librium global SFRs with only small fluctuations. However, 
the physical mechanisms are different. In the SSPSF simu- 
lations, bursts occur because percolation of star formation 
can more easily fill the simulation grid for small sizes. In the 
present simulations, the lulls of star formation in the small 
simulations occur because only a few density structures exist 
at a given time that are near the star formation threshold. 
Since star formation disrupts the local gas structure, the 
shell fragments must re-coalesce before another star forma- 
tion event can occur. Thus a "refractory time" is important 
in both models, but percolation plays no role in the present 
models. 

Figure 9 shows histograms of the star formation rate 
normalized to its mean for the last Gyr of the simulation. 
Observationally, it is equivalent to the histogram of the stel- 
lar birthrate parameter b, defined as the ratio of the cur- 
rent to past average SFR. Since the simulations studied in 
the present work are roughly statistically stationary, the 
histogram of the SFR approximates the histogram of the 
birthrate parameter b. Histograms of this birthrate parame- 
ter have been constructed by Kennicutt (1996) for a sample 
of spiral and irregular galaxies. The SFR histograms for the 
largest simulations (4 kpc) are somewhat narrower than the 
histograms constructed by Kennicutt for late type spirals or 
irregulars; however Kennicutt quotes a factor of two error in 
determining the birthrate. The long tail at large b found for 
small-size simulations is qualitatively similar to the result 
found by Kennicutt for dlrr galaxies, which are in general 
smaller than spiral galaxies. The present interpretation, in 
which the larger-b tail is due to increasing burstiness for 
smaller galaxies, is in contrast to some previous interpreta- 
tions which assume a SFR history that smoothly increases 
with time for late- type galaxies. 

The widths of these distributions are approximately a 
factor of two greater than Poission statistics would predict 
(SN = N 1 / 2 ). This disparity is not surprising since Poisson 
statistics assume independent events and spatio-temporal 
correlations of the clusters are clearly present. In fact, given 



the large coherent areas of star formation present in the 
simulations, it is perhaps surprising that the fluctuations 
are not even larger. 

4-1-2 Variation of star formation threshold 

As the star formation threshold N B f increases, the star for- 
mation rate decreases, since shells must grow to larger col- 
umn densities before forming stars, while the average num- 
ber of massive stars per cluster increases due to the assumed 
constant star formation efficiency and the larger gas mass 
available. The decrease in the star formation rate with N B { 
causes fluctuations in the star formation history to grow as 
5N cl w 27V,!/ 2 . Fi gure 10 shows the star formation history 
along with the gas velocity dispersion for three values of N s f, 
all for a domain size L = 2 kpc. For the larger values of iV s f 
(in which the star formation criterion is more stringent mak- 
ing star formation more difficult), the star formation history 
becomes bursty. 

Figure 11 shows the star formation histograms for four 
models like those presented in Fig. 9. These results suggest 
that, within the context of the present models, the observed 
distribution of present-to-average star formation rate pa- 
rameter b is sensitive to the shell column density threshold 
for star formation, besides the size effect discussed earlier. 
Thus the broad b-distributions found in dlrr galaxies could 
be due to a larger column density threshold, for example 
because of lower metal abundance. Unfortunately, this sug- 
gestion cannot be tested on a galaxy-by-galaxy basis because 
the b-distributions of the simulations refer to single models 
sampled at different times. A given small or low-metallicity 
galaxy could have a small or large value of b, depending on 
when it is observed. 

4-1.3 Time delay 

As discussed in sec. 3.2.4, increasing the time delay between 
the onset of shell instability and the onset of the cluster 
wind, Td, relative to the duration of the wind (r w = 10 7 yr), 
produces star formation structure which is spatially more 
organized into coherent "superclusters," as seen in Fig. 6. 
This enhanced spatial coherence is expected to translate into 
greater temporal coherence of the global star formation rate 
evolution. Fig. 12 shows this effect. The global SFR is shown 
as a function of time for three simulations which all have 
the standard parameters except that the time delay Td is in- 
creased from 0.5 to 2.5 to 10 Myr. Although some increased 
temporal coherence can be seen in going from Td = 0.5 to 
2.5 Myr, a large effect is seen at Td = 10 Myr, where the 
time delay equals the duration of the momentum input. Al- 
though the most evident fluctuations are seen to occur with 
timescales ~ Td, one also sees correlated behavior over all 
timescales, similar to "brown noise." The SFR does not ap- 
pear statistically stationary in the mean, although the de- 
partures from stationarity are marginal. As remarked earlier, 
larger values of the time delay are probably most appropri- 
ate if the cluster winds are driven by collective supernovae. 

It is typical of nonlinear one-zonesystems that oscilla- 
tory and even chaotic time dependence results when the de- 
lay time exceeds some other characteristic timescale of the 
system. Scalo & Struck-Marcell (1987) found, for a one-zone 
closed fluid model of star formation, that limit cycles 
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Figure 10. Evolution of global statistics for runs C with L = 2 kpc over 2 Gyr. Solid and dotted lines represent the number 
of massive stars per kpc 2 (left-hand axis) and the gas velocity dispersion (right-hand axis), rspectively. (a) t s h = 0.5 km s _1 
(N sh)21 = 3), (b) r sh = 1.5 km s" 1 (N sh)21 = 1), (c) c sh = 3 km s" 1 {N shai = 0.5). 
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Figure 11. Histograms of present-to-=past average SFR ratio for four models with different adopted column density 
threshold for star formation. The top two and bottom panels correspond to the three time histories shown in Fig. 10. 
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Figure 12. Time dependence of the global SFR (solid lines) and gas velocity dispersion (dotted lines) for simulations with 
time delay 0.5 Myr (a), 2.5 Myr (b) and 10 Myr (c). 
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(oscillations) in the system variables occurred when the time 
delay approached the characteristic timescale, which was the 
cloud collision timescale in their model. They also found 
chaotic (deterministic but unpredictable) behavior at still 
larger time delays. Unfortunately we have not yet explored 
models with Td > t w . Nevertheless, our simulations do il- 
lustrate that, even when spatial degrees of freedom are in- 
cluded, a large time delay between onset of instability and 
onset of energy injection may introduce substantial spatial 
and temporal coherence into the system, at least for the 
global scale L = 2 kpc investigated here. Whether or not 
larger time delays will result in more "bursty" behavior is 
a question left open for futue work. However we point out 
that the relative fluctuations in SFR for the Td = 10 Myr 
model are not appreciably larger than for the case with small 
Td, although the average SFR is larger, as explained in sec. 
3.2.4. 

4-1-4 Variation of average SFR with parameters 

The variation of the time-averaged SFR (per kpc 2 ) with 
model parameters is of obvious interest. Although we ran 
a large number of models, we could not explore the com- 
plete parameter space, and in most cases only varied a single 
parameter with respect to the standard model. Still, some 
interesting results can be seen from inspection of Table 1. 

First, for series A in which only the domain size was var- 
ied, the SFR (and most of the other global quantities listed 
in Table 1) are insensitive to the domain size between the 
128 2 (1 kpc 2 ) and 512 2 (4 kpc 2 ) runs. However for smaller 
sizes the SFR decreases. Notice that the "burstiness" pa- 
rameter <r(b)/b increases continuously with decreasing size, 
showing that the temporal irregularity of the SFR is truly a 
function of size, not resolution. 

The simulations of series C vary the shell internal ve- 
locity dispersion, or equivalently, the critical shell column 
density for star formation. Increasing c s h corresponds to de- 
creasing N s h,2i, so star formation takes place more easily. 
However the SFR varies roughly quadratically with the shell 
velocity dispersion for runs Cg through Cq, although the de- 
pendence is somewhat steeper for runs Ca through Cf. This 
suggests that the SFR is controlled by binary interactions 
between shells, which is plausible since, as discussed earlier, 
most of the gas accumulated by a given shell consists of 
other shells driven by different star formation sites. 

The dependence of the SFR on the time delay in series 
D follows a Td relation. Although we explained earlier why 
the SFR should increase with time delay, we can offer no 
physical explanation for the r°' 5 dependence. We note that 
if Td is inversely proportional to the square root of the den- 
sity, as should occur for gravitational instability, this gives 
a contribution to the SFR dependence ~ p _1//4 . 

The dependence of the SFR on the assumed average 
energy input per massive star (series E) is roughly linear, 
with SFR ~ E ' 8 from Table 1. 

In addition to the models discussed so far, we also ran 
one series of five 256 2 models in which the mean gas density 
was varied between values of 0.33 and 2 (the standard value 
used in all other simulations was 1.) All other parameters 
were the same as in the standard model except that the time 
delay was zero. These models are referred to as series F in 
Table 1. For the given internal shell velocity dispersion, and 



hence critical shell column density, the obvious main effect 
of increasing the average gas density is to make it easier for 
filaments to accumulate enough gas to drive them over the 
star formation threshold, increasing the star formation rate. 
Therefore this series should be similar to series C. In fact 
we do find that the SFRs of series F vary approximately 
quadratically with the mean density, again implicating shell 
interactions as the dominant process controlling star forma- 
tion. However the models of series F had zero time delay. 
From the results of series D discussed above, we expect the 
time delay to introduce an additional factor p" 1 ^ 4 to the 
SFR if the time delay varies as T d ~ 1//4 . Thus the density de- 
pendence should be somewhat shallower than p 2 . 

However there are probably "hidden variables" in the 
SFR. Scalo & Chappell (1999) present a simple model for 
the velocity dispersion scaling relation which indicates that 
the SFR should depend on, besides the parameters discussed 
above, the overall velocity dispersion of the gas c and the 
average column density through filaments \i c i. The relation 
derived, which agrees remarkably well with the present sim- 
ulations, is 

N* ~ c 3 P 2/ ' i /{E l i cl ) (16) 

Thus the density dependence of the SFR cannot be deter- 
mined observationally without explicit consideration of c 
and fi c i- The latter quantity would be especially difficult to 
determine observationally. This illustrates the danger in as- 
cribing any physical significance to an observed SFR-density 
correlation. 

4.2 Gas velocity PDFs 

Figure 13 shows the probability density functions (pdfs), es- 
timated as histograms, for three models from series C in 
which the star formation threshold is varied. The models 
are the same as for the density fields shown in Fig. 5. The 
threshold is smallest (SFR largest) for the top two panels, 
while the threshold is largest (SFR smallest) for the bottom 
two panels. The left-hand panels (a, c, e) include the entire 
density field, while the right-hand panels (b, d, f) include 
only those grid points that were identified as part of a fila- 
ment according to the filament-finding algorithm. The open 
and closed circles represent the straight and mass-weighted 
velocity fields, respectively. The histograms are presented as 
log-linear plots, so an exponential (Gaussian) pdf would ap- 
pear as a straight line (parabola). The histograms for the 
full field (left-hand panels) are based on five velocity fields 
sampled every 2 x 10 8 yr for the last Gyr of the simula- 
tions, so each composite histogram of these 256 2 runs con- 
tains 3.3 x 10 5 points. They exhibit a low-velocity core and a 
broader high-velocity component. The core appears nearly 
exponential, although a Gaussian is not excluded for the 
smallest velocities, while the broad tails have an even slower, 
possibly power-law fall-off. The broad tails reflect the veloc- 
ity distribution of the winds themselves, i.e. regions interior 
to the shells. 

Even more striking are the velocity histograms re- 
stricted to the filament ridges (right-hand panels). Since 
on average only a few percent of the lattice sites fall on 
a filament ridge, we combined filaments from 50 time steps 
equally spaced over the last half of the simulations (over 
which the global statistics appear statistically stationary; 
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Figure 13. Velocity field histograms for series C with L = 2 kpc. The shell velocity dispersions are c s h — 0.5 km s" 1 
(•Ns/,21 = 3), 1.5 km s _1 (N s f t 2i = 1), and 3.0 km s _1 (Nsf^i = 0.5) from top to bottom. The left columns are composite 
histograms of the x component of the total velocity field over five times for each simulation. The closed and open circles 
correspond to straight and mass-weighted velocities. In panels b, d, and f only those points found to lie along a filament ridge 
are included in the histograms. 
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Figure 14. Velocity histograms for 512 2 decay simulations (no star formation) with a scale-free initial energy spectrum 
E(k) oc k . Left: histogram of x-component of the total velocity field. Right: Only points found to lie along filament ridges are 
included. The parabolic reference curve corresponds to a Gaussian distribution. Comparison with Fig. 13 suggests that the 
tail excesses for the gas in filaments is not due to star formation activity. 
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Figure 15. Cloud mass functions (number of clouds per unit logarithmic mass interval) for three parameter choices in each 
of the model series' C (varying star formation threshold), D (varying time delay), and E (varying assumed energy input per 
supernova event). Each mass function was calculated by the structure tree method applied to 100 times throughout the last 
Gyr of each simulation. Solid lines at large masses have power law index -1.5. 
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see Fig. 10) to construct the histograms, resulting in about 
10 5 points per histogram. These histograms are very nearly 
exponential except for the possibly power-law flattening in 
the far tails, which is most apparent in the lowest SFR run 
(Fig. 13f). In no case is there evidence for a Gaussian com- 
ponent. 

It might be thought that the exponential pdfs are the 
result of the action of stellar winds in the simulations, which 
continually produce high- velocity filaments. In fact this is 
not the case. Instead, we propose that the exponentials are 
a consequence of momentum-conserving (but not kinetic 
energy-conserving) filament interactions. In support of this 
interpretation, we show in Fig. 14 the velocity pdf for a typi- 
cal pure decay simulation (no star formation) . The full field, 
shown on the left, is nearly Gaussian, demonstrating that 
the broad tails of the full velocity field pdfs shown in Fig. 
13 are a direct result of the star formation-driven winds. 
However the pdf of the material in filament ridges (right) 
shows a clear excess over this Gaussian for intermediate and 
large velocities. (Note that the units of the velocity axis are 
different in Fig. 14 than in Fig. 13.) The excess at large 
velocities appears exponential, even though no wind-driven 
filaments are present. The reason the exponential tails are 
not as prominent as in Fig. 13 is due to the fact that, in 
the decay simulations, the number of filaments and their ve- 
locities continually decrease with time, so the time between 
filament interactions becomes large, limiting the number of 
filament interaction which can have occurred over the dura- 
tion of the simulations. In contrast, the simulations which in- 
clude star-formation-driven winds (Fig. 13) continually sup- 
ply new high-velocity filaments, so the number of filament 
interactions which have occurred at a given (late) time is 
much greater. We therefore suggest that the exponential ve- 
locity pdfs of the filaments is a result of multiple inelastic 
shell interactions. 

Analytical models for the exponential and power law 
for tails of the velocity pdfs are given in a separate paper 
(Scalo et al. 2000), and we are mainly concerned with pre- 
senting the simulation results here. The analytical explana- 
tion for the exponential portion of the pdfs centers on the 
fact that the physical system is one in which mass and mo- 
mentum, but not energy, are the global conserved quantities. 
One way to understand this is by examining the statistics of 
a Gibbs ensemble in which there is no quadratic invariant 
(kinetic energy). Another related, but less abstract, expla- 
nation involves the statistics of momenta of interacting fila- 
ments which have undergone a number of collisions that is 
Poisson distributed. The power-law far tails for the low-SFR 
runs can be understood in terms of the predicted velocity 
pdfs of filaments or shells which have not undergone many 
interactions (see Silk 1995). At the other extreme, for a very 
large number of filament interactions, the central limit the- 
orem or a Langevin equation approach suggests a Gaussian 
pdf. Details are given in Scalo et al. (1998). 

Independent of theoretical interpretation, our simula- 
tion results are broadly consistent with the observational 
centroid velocity pdfs presented by Miesch & Scalo (1995) 
and Miesch, Scalo & Bally (1999) for regions of active inter- 
nal star formation, and with the results of optical absorption 
line and H I emission and absorption line centroid velocities 
presented in Miesch et al. (1999), both of which suggest ex- 
ponential pdfs for most regions. An exponential distribution 



of H I fluctuation velocities out of the plane of the LMC was 
found by Kim et al. (1998). Our finding of possible power-law 
contributions in the far tails of the pdf, which we attribute 
to stellar wind sources, is consistent with the analysis of the 
high- velocity optical line pdf given by Siluk & Silk (1974). 

4.3 Cloud mass spectra 

We constructed "cloud" mass spectra by employing an al- 
gorithm similar to that used in computing structure trees 
(Houlahan and Scalo 1992). The method is based on thresh- 
olding the density field at a number of intensity levels (log- 
arithmically spaced in the present analysis and given by 
Pthresh = 2' po where i is a non-negative integer) to form a 
contour image. Each closed contour is linked-up to its "par- 
ent" (defined as the next lower contour which contains the 
original contour of interest) if the parent has no other "chil- 
dren" or if the area of the original contour is larger than 
all of its siblings at that level by a predetermined factor (of 
order 2). If the region is not linked up, it is defined to be 
a cloud; otherwise it is assumed to just be part of a larger 
cloud defined at a lower threshold. 

This technique provides a sophisticated yet simple 
method of identifying clouds which avoids the problems en- 
countered by counting all closed contours over either a sin- 
gle threshold ratio or many thresholds (e.g. Stutzki et al. 
1991, Williams et al. 1994). The mass spectra constructed 
from closed density contours at a single, arbitrary contour 
level not only depend on the contour level chosen, but dis- 
agree with other techniques. In addition, methods which 
use many contour levels introduce a bias in which clouds 
are over-counted if they appear at multiple threshold val- 
ues. With the present technique, a cloud that contains no 
substructure is counted only once in the mass spectrum, in- 
dependent of the number of density contour levels chosen. 
The method is also not "fooled" by small density fluctua- 
tions that show up as small closed contours superimposed 
on larger cloud contours. Such cases may either be rejected 
as noise or counted as condensations in the larger cloud; in 
either case, the identity of the larger cloud is not affected. In 
computing the mass spectra for the simulations, we rejected 
all closed contours with areas less than five simulation cells 
from the mass spectra. 

Application of this procedure to the simulation gas dis- 
tributions produces a few hundred clouds per time step. 
Thus, to build up a statistical sample we construct a mass 
spectrum based on the density field at 100 times throughout 
the last Gyr of the simulation, giving on average 10 4 points. 
Figure 15 shows examples of the resulting composite cloud 
mass spectra for three models for three different series. (The 
mass spectra shown in the figure represent the number of 
clouds per unit logarithmic mass interval djV(m c i)/dlogm c i. 
Conversion to linear mass interval diV(m c i)/dgm c i steepens 
the power-law slope by a factor of m c \.) At large masses 
we find that all the differential mass spectra scale approxi- 
mately as diV(m c i)/dm c i oc m~j 2,5 independent of the model 
parameters. 

Comparison with observationally determined mass 
spectra is not straightforeward since the clouds derived from 
the simulation gas distribution are two-dimensional and es- 
timation of three-dimensional masses depends on assump- 
tions about the scaling between the vertical extent of the 
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clouds out of the simulation plane with their 2-D proper- 
ties. An instructive way to see that the 3D mass spectrum 
corresponding to a given 2D mass spectrum should be flat- 
ter is to assume that, given the same physical processes, 
the characteristic sizes of structures should be the same in 
3D as in 2D. More specifically, we assume that the pdf of 
structure characteristic sizes, p(r), would be the same in 3D 
as in 2D, and assume it is a power law given by p(r)~ r a , 
with a < 0. We then compare the differential mass spectrum 
f(m) — p[r(m)]dr/dm in the 2D and 3D cases. 

Inspection of simulation density images at various 
thresholds shows that most of the entities identified as 
"clouds" by the structure tree method are filamentary. We 
assume that in 3D these structures would be shells. Then, 
in 2D, the mass of a cloud is m ~ rp and the mass spectrum 
corresponding to a given size spectrum power law index a 
(assuming p = constant) is f(m) ~ m a . But in 3D we have, 
for shells, m ~ r 2 p, and obtain f(m) ~ m' Q_1 ^ 2 . Since a is 
negative, this shows that the mass spectrum in 3D will be 
flatter than the 2D simulation mass spectra by an amount 

— (a + l)/2 in the mass spectrum power law index. The flat- 
tening varies from 1/2 for a = —2 to 1 for a — —3. For 
the particular mass spectrum found in the 2D simulations, 
/(m) ~ m -2,5 , the flattening is 3/4, so the 3D mass spec- 
trum should be m _7//4 . This result is equivalent to the result 
found by simply assuming that the filaments are portions of 
shells which, in 3D, would have vertical extent comparable 
to the 2D filament length (an assumption which would not 
be justified for structures of extent larger than the vertical 
scale height). If the index of the 2D mass spectrum is 7, 
then the index in the 3D case will be (7 — l)/2, which is a 
flattening (as long as 7 < — 1). 

Additional flattening in 3D is found if we assume that 
the densities are related to size as p ~ r~ p . In that case it 
can be shown that the flattening in going from 2D to 3D is 

— (a + l)/(p 2 — 2p + 2). For any p < 2 this flattening is in 
excess of the flattening found above for the constant density 
case. For example, if p = 1/2, the 2D case with f(m) ~ 
m~ 2 ' 5 gives a = —7/4, and then the 3D mass spectrum 
should be m _3//2 . 

In reality there is not a one-to-one correspondence be- 
tween density and size in the simulations, although the cor- 
relation does exist. In the other extreme that density and 
size are independent random variables, one could in princi- 
ple use the density pdf and the assumed size pdf to calculate 
the pdf of masses m ~ rp (2D) and m ~ r 2 p (3D) using pro- 
cedures for calculating the pdf of products of random vari- 
ables. However this procedure would require the inversion of 
Mellin transforms which do not exist in closed form, and so 
we do not pursue this approach further. 

To summarize, the present 2D simulations predict dif- 
ferential mass spectra that are power laws at large masses, 
with index -2.5, independent of the simulation parameters. 
Arguments have been presented that suggest that, if the 
simulations were 3D, the mass spectra would possess large- 
mass mass spectra that are flatter by roughly unity in the 
power law index. 

A large body of observational results for molecular 
clouds find a mass spectrum slope of -1.3 to -1.9 over a large 
range of masses (see Solomon et al. 1979, Casoli et al. 1984, 
Dame et al. 1986, Blitz 1993, Brand & Wouterloot 1995, 
Kramer et al. 1998, Heyer & Tereby 1998, and references 



therein), although there are exceptions (e.g. Onishi et al. 
1996 for Taurus cores). Mass spectra with a power-law expo- 
nent of —1.5 have been derived through non-dynamical equi- 
librium coagulation models which assume the cloud velocity 
dispersions and cross sections are either mass-independent 
or have some simple form (e.g. Field & Saslaw 1965, Kwan 
1979). However the theoretical results do depend on the 
above assumptions (see Silk & Takahashi 1979, Elmegreen 
1989). For example, Elmegreen (1989) argued that if the 
cloud column densities are assumed constant in deriving the 
cross section, coagulation models should predict a steeper 
cloud mass spectrum, with n(m c i) oc m~ 2 . This result was 
used to suggest that cloud coagulation models are inconsis- 
tent with observations. (For discussions of how selection ef- 
fects could be responsible for an apparently constant column 
density, see Kegel 1989 and Scalo 1990) . N-cloud simulations 
of coalescence models (e.g. Pumphrey & Scalo 1983, Noza- 
kura 1990) can directly treat the cloud velocity distribution, 
unlike the kinetic equation approach, but still depend on 
the assumed coalescence cross section and the assumption 
of discrete independent clouds (rather than an advecting 
fluid). In contrast, the present simulations are entirely fluid- 
dynamical and involve no assumptions about cross sections 
or density-size or velocity-size relations, although they are 
limited to two dimensions. Perhaps most importantly in the 
present simulations advection can stretch or compress clouds 
into filaments, and spatial correlations between structures of 
various sizes are self-consistently accounted for. 

Our results show that the 2D mass spectra at large 
masses are approximately power laws whose slopes are insen- 
sitive to parameters, and we argue that the corresponding 
3D mass spectra are roughly consistent with observations. 
However we want to point out that our argument concerning 
the difference between 2D and 3D mass spectra is not really 
self-consistent; for example, if filaments had a vertical ex- 
tent comparable to their length, their masses and momenta 
would evolve differently, perhaps altering some essential fea- 
ture of the coagulation process. Clearly a 3D simulation with 
a careful cloud-counting algorithm will be required to esti- 
mate the 3D mass spectrum. 



5 MULTIPERSPECTIVAL INTERPRETATION 

One of the most interesting features of the present work 
is that it can be regarded as a simulation of several con- 
ceptually different models for the evolution of structures 
and star formation in galaxies, since its aspects include 
collisional coalescence of structures, gravitational instabil- 
ity in shells, propagating star formation, self-regulation, the 
dynamics of a field of advection-created discontinuities or 
shocks, and wind-driven (entirely) compressible turbulence. 
None of these aspects are more, or less, important than the 
others; they are different, incomplete, ways of viewing a 
complex phenomenon, whose complexity has already been 
severely truncated in the simulations (omission of pressure, 
magnetic fields, etc.). We feel that this "multiperspectival" 
interpretation is an important result, for it suggests that a 
particular conceptual model, while perhaps useful in describ- 
ing some aspect (s) of the observed phenomena, cannot be 
expected to provide a complete description, and will there- 
fore be "invalidated," sooner or later, by observations that 
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capture aspects excluded by the particular model. Similarly, 
when discussing conceptual models for galactic gas dynam- 
ics, a model may not be more, or less, valid than another 
model, if the two models are not meant to explain the same 
aspect. For these reasons we briefly summarize the simula- 
tions and compare with a number of previously-investigated 
conceptual models. 

In the present simulations, star formation occurs on a 
wide range of scales. First, as the wind-blown shell from a 
star formation event expands through the local gas complex 
it may trigger the birth of more stars by sweeping up the 
surrounding gas, pushing the shell column density over the 
critical value for star formation. This mode of star forma- 
tion corresponds to the scenario pictured by many propagat- 
ing star formation theories which consider a single isolated 
shell expanding into a uniform and quiescent background 
medium (e.g. Elmegreen & Lada 1977, McCray & Kafatos 
1987, Comeron & Torra 1994, Whitworth et al. 1994). How- 
ever, since the local gas complexes in the present simulations 
(and in the observed ISM) are often irregular or filamentary, 
the distribution of the triggered star formation events are 
also highly irregular in appearance, often reflecting the un- 
derlying gas distribution. This local propagation mechanism 
is found to occur in a significant fraction of star formation 
events. 

Second, star formation also proceeds through the inter- 
action of shells. When the interaction of relatively young 
shells leads to star formation, the propagation sequence is 
clear. However, star formation events involving the interac- 
tion of older shells which have undergone repeated mergings 
are the product of many initial "triggering" shells. In this 
case, star formation may also be considered as being "spon- 
taneous," although the physical process involved is the same 
as in the more local interactions. The formation of gas con- 
centrations, and thus stars, through this mechanism may be 
viewed as either taking the zero pressure limit of the hydro- 
dynamic hydrodynamic equations or including fluid advec- 
tion into a coagulation model. 

The star formation rate, controlled in part by the star 
formation threshold, is found to control the structure and 
dynamics of the gas. At low SFRs, the gas has time to coa- 
lesce into dense filaments leaving large empty regions of low- 
density gas. In addition, young shells rarely interact and the 
probability distribution of the velocity field has the power- 
law form of a set of non-interacting shells. Also, we find 
that in this case the energy spectrum is found to have the 
E(k) oc k~ 2 form predicted for a field of shocks. As the star 
formation rate increases the filaments become increasingly 
shredded and shell interactions become increasingly impor- 
tant. These interactions cause the the velocity pdf to become 
nearly Gaussian at small velocities while the tail takes on an 
exponential form. In addition, the energy spectrum slightly 
flattens as the correlation length of the shredded shells de- 
creases. 

We find that these mechanisms lead to the emergence 
of coherent and dynamic regions of star formation on spa- 
tial scales comparable to associations, complexes, and su- 
percomplexes of OB stars. The loci of these star formation 
regions move throughout the simulation, following complex 
and "unpredictable" trajectories. The spatial clustering dis- 
tribution of the young stars is scale-free, as evidenced by 



the power law correlation functions presented in Scalo & 
Chappell (2000). 

These simulations may be viewed from several perspec- 
tives. 

First, the present work may be thought of as a simula- 
tion closely related to the stellar-driven turbulence models 
of Norman & Silk (1980) and Norman & Ferrara (1996), 
who invoke severe approximations to model turbulence with 
analytic techniques. While Norman & Silk (1980) picture 
the shells as expanding into a turbulent background formed 
by a network of old momentum-conserving shell fragments, 
very similar to the present simulations, they model the tur- 
bulence through a cloud coalescence equation in which the 
basic nonlinearities governing turbulence are ignored. Nor- 
man and Ferrara (1996), on the other hand, use a transfer 
function applicable to incompressible turbulence to model 
the evolution of the energy spectrum. On the contrary, most 
of the interactions in a galactic ISM should be highly com- 
pressible. Thus, results from the present model which ex- 
plicitly model the effects of fluid advection, compressibility, 
and small-scale forcing due to stars are more general (not 
being restricted to the evolution of the correlation function) 
and correspond more closely to the supersonic conditions ex- 
pected in the ISM. However it should be emphasized that the 
results of the present simulations do correspond rather re- 
markably to the conceptual model of Norman & Silk (1980). 

Second, the present reaction-advection simulations may 
be viewed as a model of propagating star formation analo- 
gous to either the cellular automata model developed by 
Gerola & Seiden (1982, see sec. 1 for later developments) 
or the reaction-diffusion models cited in sec. 1. In our sim- 
ulations, however, the distinction between "triggered" and 
"spontaneous" star formation is not a sharp one. The pri- 
mary mode of star formation in our simulations is the colli- 
sion between dense shells or clumps. In some cases, a shell 
sweeps up and compresses enough gas to satisfy the star for- 
mation condition while the shell is still young and roughly 
circular. Such an event satisfies the traditional picture of 
propagating star formation in which a single shell expand- 
ing into a uniform ambient medium leads to the formation 
of more stars. However, the interaction of shells, after they 
have been distorted by the ambient velocity field and pos- 
sibly merged with other shells, also leads to star formation 
through the coalescence mechanism. In this case many ini- 
tial stars spread over a large region of space could be viewed 
as being equally responsible for the "triggering." 

Third, since the fluid is highly inelastic, the simulation 
may be viewed as a coalescence model which explicitly in- 
cludes the effects of fluid advection and the velocity field of 
the gas. A basic question is whether the inclusion of fluid ad- 
vection affects the form of the cloud mass spectra predicted 
by coalescence models (e.g. Field & Saslaw 1965, Kwan 1979, 
and Elmegreen 1989). Our models do not involve any as- 
sumptions concerning the form of the collision cross section 
or the "cloud" velocity distribution, both of which are com- 
pletely determined by the solutions of the reduced hydrody- 
namic equations. We find that the cloud mass spectrum at 
large masses is a power law with an index that is insensitive 
to parameters. However the precise form of the predicted 3D 
mass spectrum is uncertain. 

Fourth, the global gas structure found in our simu- 
lations, especially in the case of no star formation (de- 
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cay runs), is remarkably similar to that found in two- 
dimensional simulations of granular fluids by Goldhirsch & 
Zanetti (1993), when the assumed coefficient of restitution 
is small enough. The particles in the latter simulations are 
assumed to be disks of identical sizes which collide with geo- 
metric cross section and a given coefficient of normal restitu- 
tion. The granular system, like our galactic gas model, con- 
serves mass and momentum but not energy. Because there 
are a very large number of grains, a continuum description 
is possible (although with certain difficulties concerning the 
proper definition of thermal variables), resulting in a nonlin- 
ear advection term in the momentum equation. This simi- 
larity suggests that the overall spatial distribution is generic 
to inelastic systems that conserve mass and momentum, and 
that much of the behavior that we have found is fundamen- 
tally linked to the absence of a kinetic energy invariant. This 
idea will be used to understand the exponential velocity dis- 
tribution function of filaments in a separate paper (Scalo et 
al. 1999, in preparation). The analogy with granular fluids 
also suggests a generic relation between granular fluids and 
mass-conserving Burgers turbulence (see below). 

Fifth, we refer to this model as a "reaction-ad vection" 
model to contrast it to the reaction-diffusion systems men- 
tioned in sec. I. We suggest that a reaction-advection model 
is more appropriate to the problem of propagating star for- 
mation than reaction-diffusion models because it embraces 
more fundamental aspects of the underlying system such 
as gas conservation laws and advection. Furthermore, be- 
sides omitting important phenomena (shocks due to advec- 
tion), reaction-diffusion models have the danger of intro- 
ducing spurious phenomena (e.g. Turing patterns that de- 
pend on the dominance of diffusion). In addition, we suggest 
that generic forms of this reaction-advection model may find 
applications to problems in other fields in which advection 
plays a significant role. As far as we know, this is the first 
study of this class of systems. 

Closely related to the above view, this model can also be 
regarded as a reactive Burgers equation, since in the absence 
of stellar forcing the equations approximate Burgers' equa- 
tion (which is simply the pressureless Navier-Stokes equa- 
tion). A difference from the usual formulation of the Burgers 
problem is that we explicitly enforce mass conservation by 
solving the continuity equation as well as the (pressureless) 
momentum equation. It is interesting to note that Burg- 
ers originally studied this equation as a possible simplified 
model of Navier-Stokes turbulence, but it was later pointed 
out (Kraichman 1965) that it lacked a mechanism to de- 
stroy correlations. However, since the presence of the reac- 
tion term (star formation here) does provide this effect, we 
suggest that the reactive Burgers equation (with the inclu- 
sion of mass conservation) may be an interesting and rela- 
tively simple model of highly compressible turbulence. The 
study of "Burgers turbulence" is an extremely active field 
of research (see Bernard & Gawedzki 1998, Gotoh 1999, E 
& Vanden Eijnden 1999, Fogedby 1999, Ryan & Avellaneda 
1999 and references given there). 

6 SUMMARY 

We explored the ways in which stellar feedback and fluid 
advection can shape the morphology of the gas density 



field, the structure of the gas velocity field, and the spatio- 
temporal organization of the stars themselves. The simula- 
tions follow the evolution of a system of interacting wind- 
driven shells (which obey global mass and momentum con- 
servation laws and are subject to fluid advection). This 
allows us to investigate, in a self-consistent manner, the 
physics of wind-driven turbulence modeled by Norman & 
Silk (1980) and Norman & Ferrara (1996), propagating star 
formation models studied through the use of cellular au- 
tomata by, for example, Seiden & Gerola (1982), Comins 
(1983), Jungwiert & Palous (1994), and Pedang & Leje- 
une (1996), and diffusion approximations by, for example, 
Shore (1983), Korchagin & Ryabtsev (1992), Neukirch & 
Feitzinger (1988), and Nozakura & Ikeuchi (1984, 1988). 
An important advantage of our simulations over these mod- 
els is in our treatment of the "turbulence" through direct 
simulation of the nonlinear advection operator. One-zone 
models, by definition ignore spatial interactions altogether. 
Diffusion models include spatial couplings but in a linear 
and, we think, dangerously inappropriate manner. Norman 
& Silk (1980) model the shell-shell interactions through a 
coagulation equation which neglects, or at least highly ab- 
stracts, the nonlinear turbulent processes, while Norman & 
Ferrara (1996) use a transfer function (in Fourier space) 
for incompressible gas to model the effects of the nonlin- 
ear advection operator, and even then only at the level of 
the correlation function. Furthermore, most studies of prop- 
agating star formation ignore the gas dynamics altogether. 
In addition, our simulations, in which star formation acts 
as a "reactive" forcing term, can be compared to studies 
of Burgers' equation with external forcing. We of course do 
not claim that our simulations should serve as substitutes 
for simulations of the full hydrodynamic or MHD equations 
that include many physical processes we have excluded (e.g. 
Passot et al. 1995, who include self-gravity, the magnetic 
field, star formation feedback, heating and cooling, and ro- 
tation). However we do claim that intermediate levels of hy- 
drodynamic modeling such as that of the present paper, can 
yield much insight from an explanatory point of view, and 
also represent a more computationally efficient and struc- 
turally simpler approach than simulating the full hydrody- 
namic or MHD equations. An interesting example of this 
approach is the use of one-dimensional magnetic Burgers 
equations (Yanese 1997) to model solar flares and coronal 
heating (Galtier & Pouquet 1998). 

We now review some of the results and predictions from 
these simulations. 

(i) The spatial distribution of the gas in all the simul- 
tions evolves into a complex network of interacting partial 
shells very similar to the conceptual model envisioned by 
Norman & Silk (1980) and to the observed morphology re- 
ported by Chu & Kennicutt (1994) for 30 Dor in the LMC, 
Kim et al. (1998) for H I in the entire LMC, Mizuno et al. 
(1995) for the local Taurus complex, and Bally et al. (1999) 
for the local Circinus molecular cloud. The term "churning" 
used by Bally et al. is usefully appropriate for the present 
simulations. This dynamical network morphology seems to 
be an inescapable consequence of combined effects of stellar 
outflows, nonlinear advection, and high degree of gas com- 
pressibility. Although the inclusion of gas pressure, magnetic 
fields, and other processes will modify the results (especially 
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the thickness of the shell structures), it is hard to see how 
these processes could qualitatively alter the general mor- 
phology, which can be seen in the most physically detailed 
previous simulations (Passot et al. 1995). 

(ii) Star formation organizes into coherent structures over 
a range of spatial scales. The locus of these structures move 
throughout the grid leaving behind spatial gradients in the 
stellar ages. The sizes of these coherent structures corre- 
spond to typical sizes of OB associations, complexes, and 
supercomplexes. The correlation functions for the stars are 
scale- free power laws (Scalo & Chappell 1999). While the 
gas distributions that give rise to these structures involve a 
coagulation process, they also require feedback by star for- 
mation and the action of fluid advection. Thus, simple cloud 
coagulation models which do not incorporate such effects 
would miss this organizational mode. 

(iii) The global star formation rate per unit area is found 
to be independent of the simulation size for all but the small- 
est simulations studied. Fluctuations in the global SFR are 
found to grow in amplitude as the system size decreases. 
These fluctuations have amplitudes roughly twice that pre- 
dicted by Poisson statistics for all simulations studied and 
develop temporal correlations for simulations with spatial 
dimensions less than about a kiloparsec. Similar results were 
obtained by Gerola et al. (1980) in their Stochastic Self- 
Propagating Star Formation (SSPSF) model. They interpret 
the large fluctuations, in which long lulls occur in the global 
SFR, as as a mechanism by which dwarf galaxies can de- 
velop bursts. The present simulations, which directly model 
shell evolution and interactions and do not require the ad 
hoc assumptions of triggering and spontaneous probabilities, 
support this view. However percolation, which controlled the 
small-galaxy size bursts in the SSPSF models, plays no role 
in the present simulations. We also find that besides the 
size effect, the broad distribution of present to past average 
SFRs in late-type galaxies can be accounted for if late-type 
galaxies have smaller metal abundances, corresponding to 
a larger critical shell column density for gravitational frag- 
mentation. 

(iv) The SFR varies with the density as p 1 ' 7 if the time 
delay scales as p~ 1//2 . However we find that the SFR also 
depends on the gas velocity dispersion and on the average 
shell column density, so observed univariate correlations of 
SFR with density should be regarded with caution. 

(v) The probability distributions of filament velocities are 
strikingly exponential. The possibly power-law tails found 
at low star formation rates may reflect the power-law dis- 
tribution of a set of noninteracting shells. The origin of the 
exponential tails, however, appears to be more universal, 
since they are also found for the decay runs with no star 
formation, and we suggest that these tails are due to multi- 
ple dissipative interactions of shells. This range of forms of 
the velocity pdfs appears roughly consistent with the empir- 
ical results of Miesch et al. (1999) for local molecular clouds 
and Kim et al. (1998) for H I in the LMC. 

(vi) In the absence of star formation, we find the probabil- 
ity distribution function (pdf) of filament velocities to have 
a Gaussian core with excess, possibly exponential, tails. This 
is a new result which, to our knowledge, has not been pre- 
viously reported in the Burgers turbulence literature. The 
tail excesses are used to support our contention that the ex- 
ponential filament velocity pdfs found for simulations that 



include star formation are not primarily due to the cluster 
winds but are a result of shell interactions driven by highly 
compressible fluid advection. 

(vii) The power-law slope of the cloud mass spectrum at 
large masses is found to be insensitive to the parameters. 
The majority of the runs scale as dN(m)/dm oc m~ 2,5 in 
2D. Comparison of this result with observations is limited 
due to the two- dimensionality of the simulations. However, 
we presented arguments that the 3D mass spectrum should 
be significantly flatter, probably consistent with existing ob- 
servations. 

The central role that shell interactions play in structur- 
ing star formation in these simulations suggests that their 
detailed study could lead to a better understanding of prop- 
agating star formation. While the stability and evolution of 
single shells has been investigated through analytic meth- 
ods (Vishniac 1983, Ryu & Vishniac 1988, Vishniac 1994, 
Elmegreen 1994, 1999) and simulations (Mac Low & Nor- 
man 1993, Blondin & Marks 1996), very little work has con- 
sidered shell interactions. Stevens et al. (1992) considered 
a pair of interacting winds from a binary system, but this 
work is not directly applicable to the interaction of large, 
slower-moving shells, although it does illustrate the com- 
plexity of phenomena which may be expected. A study of 
the interaction of two shells generated by galactic explosions 
in a cosmological context (Weinberg et al. 1989) has been 
presented by Yoshioka & Ikeuchi (1990; see also Yoshioka 
& Icheuchi 1989). Although the scales of the relevant vari- 
ables (size, temperature, etc.) are much different than in the 
present simulations, the non-dimensional situation is not far 
from the conditions used here, or which might occur within 
molecular clouds. This result suggests how the present simu- 
lations may be of relevance, in a generic sense, to phenomena 
that occur in a variety of astrophysical environments from 
the cosmological to the protostellar. Analytic approaches to 
the interacting wind problem have also been developed (see 
Canto, Raga & Wilkin 1996, and references therein). We 
suggest that further detailed hydrodynamic simulations of 
the interactions between wind-blown shells in a variety of 
physical contexts is needed. In addition, the detailed effects 
of winds from young "triggered" stars embedded in an ex- 
panding shell on the shell itself would be of interest. 

Finally, we emphasize that our approach has focused on 
an intermediate-level, "reduced" description of the hydro- 
dynamics, in which several potentially-important processes 
have either been omitted (pressure, magnetic fields, heating 
and cooling, rotation) or incorporated as externally-imposed 
"rules" (gravitational instability). Our conclusions therefore 
need to be checked by three-dimensional simulations that in- 
clude these physical effects. 

This work was supported by NASA Grant NAG5-3107. 
We thank Anthony Whitworth, Enrique Vazquez- Semadeni, 
and an anonymous referee for helpful comments. 



APPENDIX A: TWO-DIMENSIONAL 
DONOR-CELL ADVECTION SCHEME 

The 2-D donor-cell advection scheme adopted for the present 
simulations was based on the methods discussed by Van Leer 
(1977) and Yanagita & Kaneko (1995). While the donor-cell 
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scheme is an Eulerian technique with a fixed grid, it is conve- 
nient to imagine that on each time step, the i th cell is trans- 
ported a distance w 1 At, where v l is the vector fluid velocity. 
Each cell donates the portion of the fluid variable which 
overlaps its downstream neighbor (see Yanagita & Kaneko 
1995 for a first-order implementation of this method). This 
method avoids the violations in certain multi-dimensional 
conservation rules known to plague schemes using the oper- 
ator splitting technique developed by Strang (1968). 

We adopt a planar approximating function for the fluid 
variable in each cell: 



A*' 3 + (x - 



x a x 



+ (y- 



\ „^7 

2/o ) a y ■ 



(Al) 



where xq and j/o are the coordinates of the cell's center and 
a];- 1 and cfy 1 are the x and y slopes. In regions of the flow 
where the velocities are oriented in the same direction, a 
cell generally receives flux from three upstream neighbors. 
The amount of donated fluid is given by the double integral 
of the approximating function over the 2-D area of overlap 
between cells. For positive velocity components, the update 
rule becomes 



A 
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Similar expressions result for velocity components with neg- 
ative or mixed signs. 

When velocities oppose one another a modification must 
be made to prevent flow crossing. We introduce this modi- 
fication by means of an example. Assume that Vi > and 
Vj+i < 0. For this example we define the gas mass which is 
to be donated from cell i to cell i + 1 to be m^. Similarly, 
mi+i is the mass donated from i + 1 to i. We stipulate that 
the fluid in cell i is successfully transferred to i + 1 only if 
it has the greater momentum, i.e. m,|«j| > m»+i|w»+i|; oth- 
erwise the opposing fluid parcel sweeps up the former and 
the mass and momentum of each are donated to the local 
cell. Thus, opposing flows are prevented from streaming past 
each other and mass and momentum are conserved. 

The slopes of the approximating functions are deter- 
mined by first order differencing (van Leer 1977) 



a: 



A i,J ~ x ) 
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(A3) 



In regions near strong discontinuities, the slopes given 
by this method can cause the approximating function to 
overshoot, leading to negative densities and/or nonphysi- 
cal oscillations. A monotonicity condition due to van Leer 
(1977) is imposed to limit the slope of the linear function, 



preventing it from exceeding the range of densities of the 
neighboring cells. In addition, if the cell is an extremum 
of the surrounding neighbors, the slope is set to zero. This 
condition is applied to the x- and y-axes independently, en- 
suring that the average value of the approximating plane 
along the cell edges is bounded by the neighbors. 



min[a\ 2\A* +1 - A % \, 2\A* -A* 1 \] 
if sgn[a l ] = sgn[A' - A' 
otherwise 



(A4) 







Near strong discontinuities, the monotonicity condition sub- 
stantially reduces the approximating slopes. A further con- 
straint is placed on the slopes to prevent any portion of 
the approximating plane from extending above or below the 
neighboring cells. 

The stability of this scheme is discussed in van Leer 
(1977) and yields the standard Courant-Friedrichs-Levy 
(CFL) condition on the velocity \v\ < 1. 

By expanding the dispersion equation to fourth order in 
kAx and second order in uAt we find the following algebraic 
dispersion relation 



(kAx) 
8At 



■[(1 



! ) - 2v(l - v)]v. 
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The steep k 4 dependence is the signature of the operator 
f4d 4 /dx 4 , where the effective diffusion coefficient U4 takes 
the form 



12 



[3(1 - V s ) - 6v(l - v)]v. 



(A6) 



Near discontinuities when the monotonicity condition 
limits the slopes a 1 , numerical diffusion is increased. In the 
limiting case where the slopes of the approximating func- 
tions are identically zero, the numerical diffusion is of the 
form V2d 2 /dx 2 , where the effective diffusion coefficient is 



"2 



IH(i-M)- 
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In each case, the numerical diffusion coefficient depends on 
the fluid velocity, reaching a maximum at v = 1/2 and going 
to zero at v = 0, 1. 

We have compared the evolution of a square wave for 
the first- and second-order schemes. We find both from an 
analysis of the simulations and a linearized analysis that 
the action of the second-order diffusion operator causes the 
length scale to evolve as I oc t 1 ^ 4 compared to the I oc t 1 / 2 
scaling for the d 2 /dx 2 operator. 

Since the model equations do not include pressure, stan- 
dard shock tube solutions cannot be used to test the code. 
Instead, we used the analytic solution of the Burgers equa- 
tion for a rarefaction shock. Details are given in Chappell 
(1997). We find that density spikes are bounded by two ac- 
cretion shocks. Between the shocks the average gas veloc- 
ity is roughly constant and equal to the the velocity of the 
shocks. The emergence of double shocks arises from diffusion 
in the continuity equation. This can be seen by including 
diffusion terms in Eq. 1 and Eq. 2 and neglecting the star 
formation forcing term. The equations become 



dp 
dt 



+ V ■ (pv) — 1/2 V p 



dpv _ _2 - 

-^j- + V • (vpv) — V2V pv 



(A8) 
(A9) 

© 1999 RAS, MNRAS 000, 



Wind-Driven Gas Networks and Star Formation 37 



where we adopt the second-order derivatives in the numer- 
ical diffusion term since we are interested in the b eha vior 
around the shocks. Substituting Eq. A.8 into Eq. A9, we 
have 
dv 
di 



+ v-Vv = v% V v + 2( V log p) (V • 



(A10) 



The first term on the right hand side represents a numer- 
ical viscosity. This term and the Lagrangian derivative on 
the left hand side of the equation constitute the Burgers 
equation. This term acts to smooth out the shock front. 
The second term on the right hand side produces the dou- 
ble shock. At the leading edge of the density enhancement, 
both the gradients of the density field and velocity field are 
negative so their product is positive and the velocity field 
is locally raised. At the trailing edge, the gradients in the 
density and velocity fields are positive and negative respec- 
tively, so their product is negative. This results in lowering 
the velocity field. Together these two effects lead to the dou- 
ble shock structure in what would be a single shock in the 
Burgers equation. We do not expect this effect to substan- 
tially influence the dynamics since the width of the density 
enhancement is found to grow only as A oc t 1 ^ 4 according to 
the fourth-order diffusion operator. Thus, the double-shock 
structure should remain as a relatively local phenomenon. 

Myhill & Boss (1993) have pointed out that this scheme 
is actually only first order accurate when the velocity field 
varies in space or time. They derive correction terms to make 
the scheme truly second order accurate. These terms effec- 
tively allow the cells to be stretched depending on the space 
and time derivatives of the local velocity field. Since we are 
interested in studying the qualitative behavior of our sys- 
tems, we do not adopt these more computationally expen- 
sive correction terms. Indeed, Yanagita & Kaneko (1995) 
have achieved remarkable success in using a truly first-order 
donor-cell scheme to study the scaling behavior and pattern 
formation in Rayleigh-Bernard convection. 



APPENDIX B: NUMERICAL DIFFUSION 

To a first approximation the diffusion along the axes is sep- 
arable and we may write the numerical diffusion coefficients 
near discontinuities as 



D 



?y = 



(Bl) 



and in smooth regions as 
1 
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The velocity dependences of these coefficients are similar, 
peaking at v = 1/2 and going to zero at v = 0, 1. 

This diffusion can be highly anisotropic. For example, a 
density spike moving along the x axis with v x = 1/2 diffuses 
along the x axis but not in the y direction. Thus, distur- 
bances moving along either axis will be stretched along their 
direction of motion. On the other hand, a particle moving 
along the diagonals with \v x \ = \v y \ will retain its shape 
since the numerical diffusion coefficients are equal. 



We have found that by adding diffusion locally only 
along the direction with the lesser numerical diffusion co- 
efficient, the anisotropy can be effectively reduced without 
substantially increasing the overall diffusivity of the integra- 
tion scheme. The amount each cell diffuses into the neigh- 
boring cells is proportional to the value of the approximat- 
ing function at the cell walls given by A\ e{t = A % — a 1 /2 
and Aright — A 1 + a*/2. Thus, the update rule for the cross 
diffusion operator may be written as 
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where e is related to the diffusion coefficient. This expression 
may be rewritten as 

1 



A 



A' + e(A l+1 - 2A l + A'- 1 ) - -e(a i 



.(B4) 



When the monotonicity condition limits the slopes to near 
zero, the last term becomes small and the diffusion takes the 
form of the standard finite difference approximation of the 
d 2 /dx 2 operator (the second term on the right hand side of 



Eq. B4). When the monotonicity condition has no effect, the 
cross diffusion takes the form 



A' 



A l + -X-A l 



-4 A 1 



■6A l + 4A* +x -A i+2 )(B5) 



where the approximating slopes (Eq. A3) have been substi- 
tuted in. This is a standard finite-difference representation 
of the d 4 /dx 4 operator. Applying Eq. B5 directly to rugged 



fluid landscapes can lead to negative densities and instabil- 
ities; however Eq. B4 subject to the monotonicity condition 



ensures a stable scheme and positive definite densities for 
e < 1/2. 

The resulting effective numerical diffusion coefficient 
along either axis may be written as 



Arff = i|Wr|[(l - M 3 ) - 2\Vr\(l - \Vr\)] 

where 



|V| 
1/2 



if \v\ < 1/2 
otherwise 



(B6) 



(B7) 



Condition B7 ensures that the diffusion coefficient is 



isotropic, depending only on the velocity amplitude. The re- 
sulting numerical diffusion increases with velocity for \v\ < 
1/2 and remains at a constant value for larger velocities. 
The transverse diffusion which must be added to the intrin- 
sic numerical diffusion of the donor-cell scheme so that the 
effective diffusion is independent of direction is given by 



D' 



D C B — D4 x 

D a g — Diy 



(B8) 



where and D^y are the numerical diffusion coefficients 
given in Eq. 22. 

The transverse diffusion terms (Eq. B8) were imple- 
mented to reduce the anisotropic numerical diffusion effects 
inherent to multidimensional donor-cell advection schemes 
of the form given in equation A3. We devised two tests to 
evaluate how well these terms preserve the shape of moving 
density disturbances. In the first test, a square density en- 
hancement was advected across the lattice at a given speed 
for a number of propagation directions. Using the unmodi- 
fied advection scheme, the disturbances were often severely 



© 1999 RAS, MNRAS 000, [i|-[m] 



38 D. Chappell and J. Scalo 



stretched along the x or y axis. We found that the inclusion 
of the transverse diffusion term preserved the shape of the 
original disturbance for all velocities examined. In the sec- 
ond test, the gas density along a circular, radially expanding 
ring was examined. In the unmodified advection scheme, the 
ring rapidly fragmented along the grid axes, while the inclu- 
sion of the transverse diffusion terms produced a substan- 
tially smoother density profile along the ring. For a detailed 
description of these and other tests see Chappell (1997). 



APPENDIX C: GRAVITATIONAL STABILITY 
OF AN EXPANDING, ACCRETING SHEET 

We consider the gravitational stability of a two-dimensional 
sheet which is accreting and expanding. This analysis fol- 
lows the work of Elmegreen (1994) and Comeron & Torra 
(1994) who were interested in the stability of a thin super- 
nova remnant or wind-blown shell expanding into a uniform 
medium. Since the present analysis slightly extends their 
work to allow independent expansion and accretion rates, 
we include some of the steps in obtaining our results. It 
is necessary to derive the more general form because the 
"environmental" density and velocity fields in our simula- 
tions quickly distort the shells, preventing the application 
of Elmegreen's instability condition which requires knowl- 
edge of the radius and age of the shell and which assumes 
a homogeneous ambient gas distribution. In this derivation, 
instabilities related to bending modes of the shell such as 
those investigated by Vishniac (1983) are not considered. 
As suggested by Elmegreen (1994), such modes may affect 
the internal velocity dispersion of the shell. 

The equations governing flows in the sheet of surface 
density /i are 



dt 



= -V ■ (jiv) + A 



dv _ _ 
H— = -VP - uv ■ Vv + ttV$ - vA 
at 

V • $ = -4?rGp. 



(CI) 

(C2) 
(C3) 



where P is the gas pressure, <& is the gravitational potential 
and A is the accretion rate onto the sheet. In addition, we 
assume an isothermal equation of state. 

We assume that the unperturbed surface density fio and 
pressure Po are uniform and that the sheet is expanding 
so that the divergence of the unperturbed velocity field is 
nonzero, i.e. V - vo / 0. Locally, however, we set vo = 0. 
The unperturbed surface density satisfies 



= A — noV ■ v . 



(C4) 



duo 
dt 

The first term on the right hand side is the rate of accre- 
tion of new material while the second term is the rate at 
which the local surface density decreases due to expansion 
of the sheet. The perturbed variables satisfy the following 
continuity and momentum equations: 



djjg 
dt 



-jUoV ■ vi - /iiV ■ v 
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fMi-gT- = -c V/ii - ^oViV ■ v - 2TTiGmfi - viA (C6) 



where we have assumed an isothermal equation of state and 
substituted V<&i = —2mGfj,i. 

Assuming the perturbations are of the form e^* - ***, the 
dispersion relation is found to be 

(u + V • vo + A/no){u3 + V • vo) = -ck 2 + 2nGfi k (C7) 

In the limit that accretion and stretching go to zero, the 
dispersion relation for the infinite sheet is recovered. The 
action of accretion and stretching do not affect the wave- 
length of the fastest growing mode since these terms do not 
contain a k dependence; they just retard the growth rate. 
Substituting the wave number of the fastest growing mode 
into the dispersion relation, we find that the condition for 
significant collapse (meaning that the growth rate exceeds a 
prescribed value co c ) becomes 



> (lj c + V ■ Vo + A//j,o)(uJc + V ■ Vo) 



(C8) 



This dispersion relation reduces to the relations found 
by Elmegreen (1994) and Comeron & Torra (1994) to within 
proportionality constants for the simpler case of a spherical 
shell expanding into a uniform medium. The effect of shell 
expansion or "stretching" (V • fib > 0) is to make collapse 
more difficult, as expected. 
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